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1 •  INTRODUCTION 


The  grid  resolution  in  any  mesoscale  meteorological  model  is  inadequate 
to  resolve  all  the  turbulent  transport  processes  likely  to  be  important  in  the 
troposphere.  Accurate  simulation  of  the  physical  mechanisms  which  are 
controlled  by  turbulence  depend  on  how  faithful  the  subgrid  flux 
parameterization  can  represent  those  turbulent  transport  processes.  In  this 
report  we  use  the  A.R.A.P.,  second-order  closure  model  of  turbulent  transport 
in  the  atmosphere  as  an  aid  in  developing  a  practical  subgrid  flux 
parameter izat ion  scheme.  The  dynamic  equations  for  the  second-order  flux 
quantities  of  interest  are  used  to  suggest  practical  approximations,  and 
relatively  high  resolution  simulation  results  from  our  existing  one  and 
two-dimensional  models  of  the  atmospheric  boundary  layer  are  used  to  test  the 
accuracy  of  tentative  schemes. 

Before  developing  a  subgrid  flux  parameterizat ion  scheme  it  is  desirable 
to  get  a  clear  concept  of  what  motions  the  parameterization  are  to  represent. 

The  parameterizat ion  should  represent  the  effect  of  the  average  of  the 
ensemble  of  motions  which  are  of  too  fine  a  scale  to  be  resolved  by  the 
mesoscale  model.  These  ensemble  average  considerations  are  discussed  in 
Appendix  A,  a  reprint  of  an  AMS  Workshop  presentation  in  October,  1983.  One 
needs  to  make  a  conscious  choice  of  the  scale  which  is  to  divide  the  fluid 
motions  into  a  resolved  mean  motion  and  an  unresolved  turbulent  motion.  The 
unresolved  ensemble  scale,  for  most  mesoscale  models,  will  include  all  of  the 
boundary  layer  scale  eddies  and  most  of  the  individual  cloud  motions.  The 
most  unique  feature  of  our  proposed  subgrid  flux  parameterizat ion  is  the 
attempt  to  incorporate  cumulus  parameter izat ion  as  an  integral  part  of  the 
turbulent  transport  model. 


We  recommend  using  a  quasi-equilibrium  approximation  to  our  second-order 
closure  transport  model.  This  approximation  combines  a  prognostic  equation 
for  the  turbulent  kinetic  energy  with  diagnostic  equations  between  local 
gradients  of  the  mean  quantities  and  the  turbulent  fluxes  and  integral  surface 
layer  relations.  Again,  the  most  attractive  feature  of  this  proposed  scheme 
is  the  incorporation  of  cumulus  parameterization  as  an  integral  part  of  the 
turbulent  transport  model. 

Results  of  six  different  tests  of  the  proposed  formulation  are  described. 
Further  testing  of  the  scheme  is  required,  particularly  for  the  cumulus 
parameterization,  but  these  preliminary  results  should  justify  a  continuation 
of  this  approach. 
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2,  PROPOSED  SCHEME 


As  discussed  in  Appendix  A,  it  is  necessary  to  ensemble  average  over 
horizontal  length  scales  which  are  at  least  4  times  the  desired  grid  length. 
Thus,  in  any  mesoscale  model  it  will  be  necessary  to  ensemble  average  over 
length  scales  which  are  much  bigger  than  the  height  of  the  planetary  boundary 
layer.  It  follows  that  vertical  gradients  of  the  mean  variables  will  be  much 
larger  than  horizontal  gradients  so  that  a  quasi  1-D  parameterization  of  the 
subgrid  turbulent  fluxes  is  appropriate;  i.e.,  although  the  variables  are  a 
function  of  time  and  all  3  space  coordinates,  the  subgrid  fluxes  will  be 
primarily  driven  by  the  vertical  gradients.  After  reducing  the  problem  to  a 
quasi  1-D  problem,  it  is  still  necessary  to  decide  on  the  level  of  the 
turbulent  closure  to  be  used.  Although  full  Reynolds  stress  closure  may 
eventually  prove  desirable,  a  reasonable  compromise  between  the  desire  to 
provide  as  much  physics  in  the  turbulent  transport  as  possible  without  unduly 
increasing  the  numerical  complexity,  .  is  to  use  a  quasi-equilibrium 
approximation  for  the  turbulent  flux  equations.  This  approximation  combines  a 
prognostic  equation  for  the  turbulent  kinetic  energy  with  diagnostic  relations 
between  local  gradients  of  the  mean  quantities  and  the  turbulent  fluxes. 

The  equation  for  the  turbulent  kinetic  energy  may  be  modeled  as 
( Lewellen ,  1 981 ) : 


D(q2/2) 

Dt 


u'w'  —  -  v'w’ 
Sz 


(2.1) 


The  last  2  modeled  terms  representing  diffusion  by  turbulent  transport  and 
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turbulent  dissipation  are  not  specified  until  the  turbulent  length  scaie  A  is 
specified.  Again  consistent  with  our  desire  to  minimize  increased  compxexity 
we  choose  the  following,  relatively  simple  algebraic  specification: 

A  -  min  (0.65  z,  A*,  q/2N)  (2-2) 

unless 


dA 

dz 


>  0.65,  in  which  case 


I  dA 
I  dz 


0.65. 


(2.3) 


z  is  height  above  the  surface,  q^  ■  uiui  twice  the  turbulent  kinetic 
energy,  and  N2  -(g/T0)  dT/dz  is  the  Brunt-Vaisala  frequency  if  positive, 
otherwise  this  term  is  neglected  in  the  minimization.  A*  is  a  measure  of  the 
scale  of  the  turbulent  region  and  is  set  equal  to  O.22Z4,  where  Z4  is  the 
height  at  which  q2  becomes  one  quarter  of  the  surface  value,  i.e.,  the  outer 
edge  of  the  boundary  layer.  The  constant  0.22  was  chosen  to  obtain  good 
agreement  with  experimental  observations  of  w'w'.  The  second  restriction  on 
the  slope  of  A  is  used  to  prevent  A  from  decreasing  too  sharply  in  any 
inversion  region,  since  it  would  be  physically  unrealistic  for  the  lengthscale 
to  change  too  abruptly.  The  scale  may  be  determined  by  sweeping  upward  from 
the  surface  applying  the  minimum  criterion  until  the  slope  bound  is  exceeded; 
then  the  scale  is  set  using  the  scale  at  the  point  below,  together  with  the 
slope  criterion.  This  simple  specification  of  A  is  probably  one  of  the  major 
limiting  factors  on  our  proposed  parameterization,  but  any  of  the  available 
more  complicated  formulations  do  not  appear  warranted  for  this  level  of 
approximation  at  the  present  time. 


The  diagnostic  relations  for  the  vertical  turbulent  fluxes  are  obtained 
by  assuming  local  eQuilibrium  in  all  of  the  equations  involving  second  orde 
correlations  except  the  horizontal  velocity  fluctuations.  This  yields. 


u 1  w r 


QASn 


3u 

3z 


(2.4) 


2-2 


(2.5) 


»  U  t 


v  1  w 


w’e*  - 


36 

-  qASHiI 


(2.6) 


3h 


Wh'  -  -  qASh  ^ 

where  Sm  and  SH  are  the  following  stability  dependent  coefficients: 

n  ( 1 ~2b) _ 


(2.7) 


H  3CA  +  (2  +  1/bs)R] 


(2.8) 


A  ♦  d/bs  -  1/A)R  (2.9) 

m  H  1  ♦  r/a 

where  R  -  (gA2/q2T0 ) 3T/3z,  and  A,b,  and  s  are  our  standard  model  coefficients 
(0.75,  0.125,  1.8). 

Equations  2.8  and  2.9  are  somewhat  simpler  than  the  expressions  used  in 
Mellor  and  Yamada's  (1974)  level  2-1/2  model  because  we  have  assumed  the 
vertical  velocity  variance  equation  to  be  in  local  equilibrium  while  their 
relations  assumed  1/3rd  of  the  local  q  imbalance  in  Equation  2.1  to  be 
apportioned  to  the  wfwf  balance.  Which  is  more  physical  is  arguable  since  an 
imbalance  under  atmospheric  boundary  layer  conditions  appears  more  likely  to 
persist  in  the  horizontal  velocity  variance.  Our  principal  reason  for 
choosing  local  balance  in  wf  wf  to  close  the  equations  is  because  this  leaves 
Sm  and  independent  of  the  local  velocity  gradient. 

This  quasi-equilibrium  parameterization  of  the  subgrid-scale  fluxes 
provides  a  robust  formulation  except  under  conditions  which  may  lead  to 

extreme  values  of  R.  Under  unstable  conditions  SH  has  a  singularity  near 
R  =  -0,12.  Generally,  we  do  not  expect  this  limit  to  be  reached  since  as  SH 
becomes  large  the  heat  flux  production  term  in  Equation  2.1  increases  q2  and 
thus  reduces  R.  For  added  robustness,  S ^  is  not  allowed  to  exceed  2,  in  which 
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case  Sm  is  also  set  equal  to  2.  This  limit  is  set  high  enough  it  is  rarely 
invoked,  but  low  enough  to  not  cause  problems  when  it  is.  Also,  R  is  not 
allowed  to  exceed  1,  which  gives  a  lower  limit  of  0.035  for  SH;  note  that 
this  is  only  a  restriction  in  the  region  where  the  slope  criterion  governs  the 
choice  of  A,  since  otherwise 


A  <  q/2N,  i.e.,  A 


2  _£  H 
„  3z 


q2/4 


<2.10) 
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3.  SURFACE  BOUNDARY  CONDITION  PARAMETERIZATION 


Appropriate  surface  boundary  conditions  for  any  numerical  meteorological 
simulation  depend  critically  on  the  vertical  grid  resolution  next  to  the 
surface.  We  wish  to  provide  a  formulation  of  the  surface  transport 
coefficients  which  does  not  require  the  detailed  resolution  of  a  part  of  the 
surface  layer  flow.  We  would  like  our  formulation  to  be  completely  compatible 
with  the  standard  surface  layer  relationships  but  yet  sufficiently  robust  that 
it  can  give  approximately  valid  results  when  the  minimal  resolution  is  so 
coarse  as  to  include  the  total  boundary  layer  within  the  model  layer  adjacent 
to  the  surface. 

We  choose  to  define  the  surface  transport  coefficients  by: 


Cf  =  -  u'w' 
ru 

0  1 

f  qu1 

(3.1 ) 

* 

1  > 

1 

II 

> 

o 

'oAvi 

(3-2) 

C0  =  -  w * Q' 

'°  i 

o 

o 

CT 

(3-3) 

when  the  subscript  o  represents  the  surface  value  at  zQ  and  the  subscript  1 
represents  the  average  of  the  quantity  between  the  surface  and  midway  between 
the  first  two  grid  points.  We  designate  this  height  as  za. 

We  choose  to  involve  q,  in  the  definitions  of  Equations  3-1 "3-3.  rather 
than  use  either  u,  -or  u*  in  the  place  of  this  characteristic  velocity  as  is 
most  often  done.  This  removes  some  of  the  variation  from  the  transport 
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coefficients  without  adding  complications  since  the  surface  value  of  q  already 
is  needed  as  a  boundary  condition  on  Equation  2.1.  With  this  definition  the 
coefficients  are  not  forced  to  become  very  large  as  free  convection  conditions 
are  approached.  The  boundary  conditions  are  now  completed  by  providing  a 
specification  for  the  surface  transport  coefficients  and  the  turbulent 
velocity  fluctuation  q.  An  expression  for  q^  is  derived  by  integrating  the 
turbulent  kinetic  ene'rgy  equation  between  z0  and  za. 
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For  purposes  of  our  surface  parameterization  we  will  ignore  the  LHo  of 
Equation  3.4  since  the  time  scale  of  the  turbulence  in  the  surface  layer 
should  be  short  in  comparison  to  the  time  scale  of  the  mesoscaue  model 
simulation.  The  turbulent  scale  A  appearing  in  the  last  term  can  be  taken  as 
proportional  to  z,  i.e.  az,  up  to  some  height  which  depends  on  either 
stability  or  the  boundary  layer  thickness.  A  reasonable  representation  of 
Equation  3 • 4  then  should  be  to  take 
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where 
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and  zs  is  a  height  based  on  the  maximum  value  of  A  which  is  yet  to  be  defined. 
After  some  experimentation  with  this  expression  we  found  it  desirable  a.so  to 
limit  the  heat  flux  term  to  damp  out  no  more  than  20$  of  the  shear  flow  term 
under  stable  conditions.  This  corresponds  to  imposing  an  upper  bound  on  the 
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flux  Richardson  number. 


Approximations  for  the  transport  coefficients  may  be  obtained  from 
integrals  of  the  Reynolds  stress  and  heat  flux  equations  between  z0  and  za. 
We  begin  with  the  Reynolds  stress  equations  as  modeled  by  our  second-order 
closure  approach( Lewellen ,  1981)  and  integrate  them  between  z0  and  za.  Under 

neutral,  homogeneous,  quasi-steady  conditions,  these  equations  may  be  written 
as: 
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When  A  is  assumed' equal  to  az,  q  and  w'w'  are  taken  as  constant,  and  u'w', 
v'w',  and  Q'w'  are  taken  as  linear,  these  equations  reduce  to 
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When  the  turbulent  fluxes  at  za  are  specified  by  the  subgrid  flux 

parameterization  between  points  1  and  2,  the  mean  value  of  u,v,  and  0  are 

-  2 

averaged  between  their  respective  values  at  points  1  and  2,  and  w’w'/q  is 
approximated  as  equal  to  1/4;  Equations  3 - T  0—3  - 1 2  provide  approximate 
expressions  for  the  determination  of  the  surface  transport  coefficients 
defined  in  Equations  3-1 “3 -  3 - 

If  za  falls  within  the  surface  layer  where  u'w'a  -  u'w'0,  u  -  In  z/z0, 
etc.,  these  last  3  equations  reduce  to 
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Where  we  have  given  a  wfwf/q  its  neutral  value  of  0.16* 


(3.1*0 


There  are  3  types  of  modifications  which  must  be  considered  for  Equations 

3.13  and  3.14.  There  are  the  stability  modifications  which  can  be  made  within 

the  framework  of  the  constant  flux  surface  layer.  There  are  the  effects  of 

pressure  gradients  and  heat  sinks  which  force  the  turbulent  transport  fluxes 

to  depart  from  their  surface  values.  Finally,  there  are  modif icat ions 

required  when  z  extends  above  the  height  to  which  the  turbulent  scale  can  be 
a 

taken  as  proportional  to  z.  Let  us  consider  stability  modifications  first. 

The  surface  layer  relationships  are  familiarly  given  in  terms  of  z/L. 
i.e.,  for  stable  conditions 
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Equations  3.15  and  3.16  may  be  integrated  and  inverted  to  give  expressions  for 
the  transport  coefficients  as  a  function  of  za/L.  In  order  to  Improve  the 
robustness  of  the  formulation  we  choose  to  approximate  the  stability  factor  in 
terms  of  the  layer  bulk  Richardson  number. 
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When  Equations  3.15  and  3-16  are  integrated  and  order  zQ/za  neglected  it  is 
possible  to  write 
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An  approximation  for  the  stable  modification  which  agrees  very  well  with 
the  standard  surface  layer  relation  is 
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The  2  expressions  for  cf/cfN  in  Equation  3-20  are  compared  in  Figure  3-1-  The 
two  forms  are  almost  indistinguishable.  The  compatible  expressions  for  Cq  are 
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These  are  compared  in  Figure  3-2.  We  expect  the  Rig  formulation  to  be 
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FIGURE  3.1:  COMPARISON  OF  THE  TWO  EXPRESSIONS  FOR  Cf/Cf  IN  EQUATION  3.20  AS  A 

n 

FUNCTION  OF  BULK  RICHARDSON  NUMBER. 
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FIGURE  3.2: 


COMPARISON  OF  THE  TWO  EXPRESSIONS  FOR 
FUNCTION  OF  BULK  RICHARDSON  NUMBER. 
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approximately  valid  even  when  a  coarse  grid  is  used  and  the  model  surface 
layer  extends  well  beyond  the  actual  surface  layer.  However,  in  this  case,  za 
can  no  longer  be  determined  solely  by  the  grid  but  also  must  be  a  function  of 
the  boundary  layer  thickness  which  will  impose  an  upper  bound  on  Rig  to 
prevent  Equations  3.20  and  3-21  from  becoming  negative.  This  will  be 
discussed  later. 

For  unstable  flow  we  propose  to  use  essentially  the  same  formulation.  We 
do  not  anticipate  appreciable  negative  values  of  Rie  being  maintained  for 
significant  times  since  the  strong  instability  should  generate  Increased 
turbulence  which  tends  to  wipe  out  the  unstable  gradients.  This  change  in 
transport  is  already  accounted  for  by  the  presence  of  q  in  Equations  3.1  and 
3.2.  For  robustness,  we  will  not  let  -Rig  exceed  1/4  which  is  about  as  far  as 
the  linear  form  of  Equation  3-17  should  be  expected  to  hold,  and  to  add  q2, 
twice  the  turbulent  kinetic  energy,  to  the  denominator  of  Equation  3.17  so 
that  it  remains  limited  even  under  free  convection  conditions.  The  specific 
limit  of  free  convection  is  considered  in  Appendix  B. 

Let  us  now  look  at  the  effects  of  pressure  gradients  and  heat  sinks  which 
force  the  turbulent  transport  fluxes  to  depart  from  their  surface  values. 
This  extension  of  the  neutral  surface  layer  relations  beyond  the  region  of 
constant  flux  follows  the  analysis  given  by  Lewellen  (1977).  At  the  surface 
where  the  advective  terms  are  zero  the  mean  momentum  and  mean  energy  equations 
may  be  approximated  as 


3  u’w’  _  F  (3-22) 

“aT"  1 


l  Q,w>  -  q  (3-23) 

9z 

where  F-  is  the  pressure  gradient  vector  and  Q  represents  any  distribution  of 
heat  sources  (or  sinks).  When  Equations  3-22  and  3-23  are  integrated  from  z0 
to  z  they  allow  the  turbulent  fluxes  at  za  to  be  related  to  those  at  the 
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surface.  Thus 
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These  last  3  equations  can  be  used  in  Equations  3*  1 0-3*1 2  to  obtain  the 
influence  of  F  and  Q  on  the  surface  transfer  coefficients.  We  use  these 
expressions  to  find  approximate  expressions  for  and  Cq^,  namely 


Equations  3*27  and  3-28  are  appropriate  as  long  as  za  is-  less  than  the 
height  at  which  the  pressure  gradient  term  balances  the  surface  shear  stress. 
If  these  balance  heights  are  defined  as  z*,  then  they  are  given  from  Equations 
3. 24  to  3.26  as: 

*  ,  /  1  3P  | 

*u-  /pi tj  <3.29) 


z 


# 

V 


/ 


1  JE_ 

p  3y0 


(3.30) 


* 


ze  “ 


(3-3D 


3-9 


Now  we  are  ready  to  face  the  problem  of  extrapolating  the  formulations  to 
the  outer  regions  of  the  boundary  layer  where  A  no  longer  increases  with  z  and 
where  the  turbulence  can  not  be  assumed  equal  to  its  surface  value.  These  two 
effects  tend  to  partially  cancel  each  other.  In  order  to  account  for  these 
effects  in  any  definitive  way  it  would  be  necessary  to  estimate  the  thickness 
of  the  boundary  layer.  We  have  experimented  with  integrals  of  the  mean 
momentum  equations,  the  temperature  equation  and  the  turbulent  kinetic  energy 
to  provide  estimates  of  the  boundary  layer  thickness  when  the  vertical 
resolution  of  the  model  is  too  coarse  to  provide  a  valid  estimate.  This 
formulation  is  presented  in  Equations  3*35  to  3*51,  but  the  tests  results  to 
be  presented  in  the  next  section  suggests  that  a  fairly  consistent  formulation 
is  possible  without  adding  the  complexity  of  prognostic  equations  for  these 
boundary  layer  thicknesses.  This  simpler  formulation  is  obtained  by  limiting 

|Rib|  £  0.2  (3*32) 


and  under  stable  conditions  to  set 


6  *  0.2 


(3*33) 


For  purposes  of  estimating  zs  in  Equation  3*6>  the  maximum  A  may  be  estimated 
by  setting  the  stable  boundary  layer  thickness  as  equal  to  5L 
(Deardorff,  1972),  and  Amax  «  0.2L  (Lewellen  &,Teske,  1973).  This  yields 


z  c  =  0.012 


(3-34) 


Equations  3.20  and  3.21  are  extrapolated  beyond  the  surface  layer  by 
restricting  the  numerator  to  be  limited  as  if  Rig  -  0.2,  but  not  limiting  Rig 
in  the  denominator.  On  the  unstable  side  Rig  is  everywhere  limited  to 

<  0.2. 
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The  effect  of  the  pressure  gradient  is  extended  by  replacing  the  surface 
pressure  gradient  with  f  (vg  -  v,)  and  f  (u1  -  ug)  in  Cfu  and  Cf^  respectively. 
We  have  not  attempted  to  include  surface  layer  radiation.  For  robustness  we 
have  also  included  the  restriction  that  the  angle  of  the  surface  shear  stress 
can  not  depart  from  arctan  v^  /  ui  by  more  than  45°. 

This  completes  the  formulation  required  for  the  tests  in  the  next 
section. 


On  the  unstable  side  the  tests  provided  in  the  next  section  generally 
have  sufficient  grid  points  to  provide  a  rough  estimate  of  5.  If  the 
formulation  is  to  be  used  with  a  much  coarser  resolution  then  even  the 
unstable  layer  is  likely  to  appear  as  a  stable  layer.  For  this  coarser  grid 
resolution  it  would  be  desirable  to  include  an  integral  equation  for  6.  A 
number  of  prognostic  equations  for  6  have  been  published  e.g., 
Deardorff,  1972,  Randall,  1979.  A  consistent  formulation  of  6  may  be  derived 
from  integrals  of  the  mean  momentum,  temperature,  and  turbulent  kinetic  energy 
equations.  The  momentum  and  temperature  equations  may  be  written  in  the 
following  defect  form 


—  u'w'  -  f(v 
3z  ' 


(3-35) 


3 

3t 


(v.  -  v) 


-i  v’w*  +  f/U 
3z  V 


)* 


_3  v 

at 


(3-36) 


3 

3t 


$  +  —  0 
at  0 


(3-37) 


where  by  definition  the  subscript  ®  represents  the  value  of  the  quantity  in 
the  inviscid  region  above  the  boundary  layer,  fVg  and  fUg  represent  the  two 
components  of  the  pressure  gradient,  and  <2  is  any  thermal  energy  source  term. 
We  can  replace  V  and  U  in  Equations  3-35  and  3-36  with  and  U^,  after 
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after  integrating  across  the  boundary  layer  Equations  3 • 35 ~3 • 3T  may  be 
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The  turbulent  kinetic  energy  may  be  written  as: 


When  we  define 


(3-49) 


and  assume  the  turbulent  fluxes  vary  linearly  across  the  boundary  layer 
Equation  3.48  can  be  rewritten  as: 
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A  reasonable  form  for  the  dissipation  appears  to  be 
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where  A  and  B  are  constants  to  be  obtained  by  fitting  high  resolution  runs. 


During  high  resolution  runs  none  of  these  last 
be  required  but  for  the  general  case  we  need 
equations  for  6U,  6yf  6g,  and  q^  and  provide  a 
determining  6. 


equations  3 - 40—3 * 51  should 
to  carry  these  prognostic 
consistent  algorithm  for 
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We  would  like  5  to  be  a  valid  measure  of  the  height  to  which  the 
turbulence  in  the  boundary  layer  extends.  Under  neutral  conditions  this  can 
be  provided  by  a  simple  proportionality  constant,  i.e., 


(3-52) 


However,  when  strong  thermal  gradients  are  present  the  spread  of  the 
turbulence  is  more  likely  to  be  controlled  by  the  thermal  thickness  than  it  is 
by  the  velocity  displacement  thickness.  A  practical  scheme  for  accomplishing 
this  transition  will  require  some  numerical  experimentation.  The  unstable 
transition  is  the  best  understood.  The  inversion  at  the  top  of  the  boundary 
layer  takes  over  the  governing  role,  with  the  rate  of  entrainment  across  this 
inversion  determining  the  growth  rate  for  6.  Under  these  conditions,  when  9 
is  nearly  uniform  across  the  boundary  layer,  Equation  3-45  reduces  to 


(3.53) 


with  A©^  equal  to  the  temperature  change  across  the  inversion.  This 
temperature  jump  may  be  specified  or  may  be  estimated  by  setting  the  Ri  across 
the  inversion  equal  to  a  critical  value  arid  assuming  its  thickness  to  be  a  set 
fraction  of  the  boundary  layer  thickness. 

If 
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is  fixed,  then 
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(3.55) 


To  begin  with  it  is  safer  to  specify  A©j_ ,  but  Equation  3-55  should  provide  a 
reasonable  rough  estimate  with  Ri^/Y  *  1.2.  With  A0j_  fixed  Equation  3.53  can 
be  used  to  eliminate  5q  and  the  temperature  equation,  Equation  3- ^2,  becomes  a 
prognostic  equation  for  6.  A  number  of  formulations  of  this  equation  under 
these  conditions  have  been  given  in  the  literature. 


The  thermal  control  of  the  boundary  layer  under  stable  conditions  is  more 
subtle.  Under  stable  conditions  we  expect  the  velocity  gradient  and  the 
temperature  gradient  to  be  constant  over  the  majority  of  the  boundary  layer 
between  the  surface  and  the  top  of  the  boundary  layer.  The  essential 
dynamical  control  should  be  obtained  by  limiting  the  Ri  based  on  these 
gradients  to  a  critical  value. 
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Under  stable  conditions  we  expect  the  turbulent  scale  to  be  limited  by  the 
ratio  of  the  turbulent  energy  to  Brunt-Vaisala  frequency,  i.e. 
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In  this  case  Ri0  can  be  computed  from  Equation  3-18  with 
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(3.59) 


Thus  6  is  computed  as  the  minimum  of  Equation  3-52  .or  Equation  3*57  and  Rig  is 
the  minimum  of  Equation  3.17  with  6  replacing  za  or  Equation  3-18  with 

Equation  3-59. 


In  summary,  the  surface  boundary  conditions  to  be  used  with  the 
quas i-equil ibri urn  formulation  of  section  II  in  the  tests  to  be  discussed  in 
Section  V  are: 
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The  system  is  closed  by 
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unless  a  separate  equation  is  required  for  6. 
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4.  RESULTS  FOR  SOME  TEST  PROBLEMS 


The  best  test  of  any  sub-grid  flux  parameterization  is  in  how  well  it 
allows  the  '.model  within  which  it  is  used  to  agree  with  real  world  data. 
However,  it  is  tough  to  carry  out  conclusive  tests  of  this  type  for  two 
reasons.  First,  it  is  difficult  to  separate  out  the  effects  of  the  flux 
parameterization  from  the  effects  of  the  remaining  model  characteristics  and 
second,  it  is  seldom  possible  to  find  sufficiently  complete  data  that  all  the 
uncertainty  involved  in  input  conditions  can  be  removed.  We  have  chosen  an 
easier  type  of  test,  which  involves  comparing  the  results  of  coarse  grid, 
highly  parameterized  simulations,,  with  high  resolution,  less  parameterized 
simulations. 

We  use  the  results  of  high  resolution  simulations  using  our  full 
second-order  closure  model  of  turbulent  transport  as  a  standard,  and  wish  to 
show  the  extent  to  which  our  low  resolution  parameterized  model  proposed  in 
Sections  2  and  3  will  agree  with  these  high  resolution  results.. 

All  the  results  presented  in  this  section  come  from  1-D  simulations  where 
only  vertical  gradients  are  considered.  These  are  appropriate  tests  for 
sub-grid  parameterization  for  a  model  like  NORAPS  where  the  horizontal  grid 
resolution  is  likely  to  be  approximately  100km.  However,  to  support  the 
results  of  our  1-D,  high  resolution,  full  second-order  closure  model  we  have 
also  conducted  some  tests  using  a  2-D  version  of  our  full  mcdex.  The 
advantage  of  the  2-D  simulations  is  that  they  are  made  under  conditions  when 
we  expect  much  of  the  turbulent  transport  to  be  carried  by  relatively  large 
2-D  eddies,  which  can  be  resolved  by  the  grid  system  used.  Thus,  the  results 
show  little  sensitivity  to  the  closure  assumptions.  These  2-D,  large  eddy 


simulations  then  provide  direct  tests  of  the  sub-grid,  flux  parameterization 
used  in  the  high  resolution  simulations.  The  results  presented  in  Appendix  C 
show  that  the  1-D,  high  resolution  results  are  quite  close  to  the  2-D,  LES 
results,  even  under  the  unstable  conditions  where  simple,  second-order  closure 
results  have  been  reported  to  fail  to  provide  satisfactory  entrainment  results 
(Zeman  and  Lumley,  1976),  Wyngaard,  1981).  These  2-D  results  are  somewhat  of 
a  digression  from  our  main  task  of  providing  sub-grid  parameterization 
appropriate  for  a  mesoscale  model,  but  we  felt  it  was  necessary  to  provide 
clear  support  for  the  accuracy  of  the  1-D  full  closure  equations  before  using 
these  results  as  a  standard  for  our  parameterization  tests. 

We  have  chosen  3  quite  different  meteorological  regimes  to  test  our 
proposed  parameterization.  The  first  is  the  unstable,  relatively  well  mixed 
planetary  boundary  layer  under  conditions  quite  similar  to  those  used  in  the 
2-D  simulations  of  Appendix  C.  The  second  regime  is  that  where  a  constant 
cooling  rate  is  applied  which  allows  the  PBL  to  approach  an  equilibrium  stable 
layer.  The  third  test  is  where  we  allow  the  surface  heat  flux  to  transition 
from  unstable  to  stable  and  back  to  unstable  on  a  cycle  which  simulates  a 
typical  diurnal  variation.  The  full  closure  model  was  run  with  a  sufficiently 
fine  grid  to  make  the  solution  independent  of  the  grid.  The  Q.E.  model  was 
run  with  4  different  grid  resolutions  as  given  in  Table  1. 

Comparisons  between  the  Q.E.  results  and  the  full  closure  results  for 
these  three  tests  are  given  in  Figures  4.1  to  4.17.  When  the  great 
simplification  in  the  model  is  considered,  the  Q.E.  results  appear  quite 
acceptable. 

Results  for  an  unstable  PBL  growing  into  a  stable  region  similar  to 
conditions  of  Appendix  C  are  given  in  Figures  4.1  to  4.2.  The  Q.E.  result 
shows  a  negative  temperature  gradient  across  the  middle  of  the  boundary  layer 
which  is  not  present  in  the  results  of  the  full  closure  model.  This  gradient 
is  required  to  maintain  the  positive  heat  flux  across  the  bulk  of  the  layer. 
The  diffusion  terms  in  the  full  closure  model  allow  the  heat  flow  with  little 
or  no  gradient  across  much  of  the  layer.  In  spite  of  this  difference  the 
results  are  reasonably  close.  Results  for  the  free  convection  simulation  of 


Figures.  4.4  and  4.5  are  quite  similar. 

Results  for  the  constant  cooling  rate  case  are  shown  in  Figures  4.3 

-4  -1 

to  4.7.  A  nominal  neutral  PBL  with  Ug  *  10  m/sec,  zQ  *  0.01  m,  f  *  1 0  -  — 

is  simulated  for  12  hours  and  then  a  surface  cooling  rate  of  1.8  °C/hour  is 
applied.  The  near  surface  Az  for  each  of  these  is  indicated  on  the  figure. 
Figure  4.3  shows  that  the  Q.E.  model  tracks  the  surface  friction  velocity 
very  well  for  all  of  the  grids  tested.  The  surface  heat  flux  comparison  dues 
not  show  quite  as  much  agreement.  The  larger  the  grid  size,  the  more 
disagreement  over  the  first  3  hours.  However,  this  improvement  with  grid 
resolution  is  not  consistent  in  the  asymptotic,  late/time  results. 

The  vertical  profile  of  u,  v,  and  T  are  given  in  Figures  4.5  and  4.6  at 
one  particular  time  6  hours  after  the  cooling  began.  Even  the  coarse  grid 
Q.E.  results  are  quite  reasonable.  The  maximum  turbulent  velocity  is  shown 
in  Figure  4.7.  The  full  closure  result  is  significantly  higher  than  the  Q.E. 
results  principally  because  the  maximum  q  is  occurring  at  the  surface  and  the 
Q.E.  models  use  a  turbulent  velocity  averaged  over  a  significant  fraction  of 
the  boundary  layer. 

Results  for  the  diurnal  day  variation  are  shown  in  'Figures  4.8  to  4.17. 
The  largest  discrepancies  occur  in  the  morning  transition  when  there  is  a  thin 
unstable  layer  eroding  the  surface  inversion  which  has  built  up  over  the 
night.  This  causes  a  big  difference  between  the  course  resolution  and  fine 
resolution  Q.E.  results  over  approximately  6  hours. 

Other  comparisons  are  shown  in  Section  V  and  in  Appendices  C  and  D. 
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TABLE  1 


Computational  Grids  -  Q.E.  Model 


Constant  Cooling  -  Uniform  grid  spacing 


Grid 

Az 

No.  Points 

z  (1) 

A 

140 

7 

70 

B 

100 

10 

50 

C 

70 

14 

35 

D 

30 

31 

15 

Diurnal 

,  Day  and  Unstable  Boundary  Layer 

Uniform  grid  spacing 

Grid 

Az 

No.  Points 

z  (1) 

C 

140 

7 

70 

Nonuniform  grid  spacing 

Grid 

Az  (1) 

No.  Points 

z  Cl) 

A 

140 

31 

70 

B 

60 

39 

30 

Unstable  Boundary  Layer 


FIGURE  4  1-  COMPARISON  OF  THE  POTENTIAL  TEMPERATURE  DISTRIBUTION  AS  PREDICTED 
BY  THE  Q.E.  MODELS  OF  TABLE  1  WITH  THAT  PREDICTED  BY  THE  FULL 
CLOSURE  MODEL  FOR  A  QUASI-STEADY  UNSTABLE  FLOW. 
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Unstable  Boundary  Layer 


FIGURE  4.2:  COMPARISON  OF  THE  HEAT  FLUX  DISTRIBUTIONS  WHICH  GO  WITH  THE 
POTENTIAL  TEMPERATURE  DISTRIBUTIONS  OF  FIGURE  4.1 
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Constant  Cooling  Rate 


FIGURE  4.3:  FRICTION  VELOCITY  AS  A  FUNCTION  OF  TIME  AFTER  A  CONSTANT  COOLING 
RATE  OF  1.8°c/H0UR  IS  APPLIED  TO  THE  SURFACE. 
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Surface  Heat  Flux 


Constant  Cooling  Rate 


FIGURE  4.4:  SURFACE  HEAT  FLUX  AS  A  FUNCTION  OF  TIME  AFTER  A  CONSTANT  COOLING 
RATE  OF  I.8°c/H0UR  IS  APPLIED  TO  THE  SURFACE. 
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Constant  Cooling  Rate 


r 


T 


T 


T - 1 - 1 - 1 - J - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1  I  I  *  *  <  1  I  r 


A  -  140  A  — 

A  -  100  — - B  — 

A  -  70  C  — 

A  -  30  0  — 


Full  Closure  Model 


,  .  f  .  i  t  t  i  t  i  »  i  i  < _ i _ i _ i — l — i — i — i — i — i — i — » — i — i — t- 

12  18  24  30 

Time  (  Hours  ) 

MAXIMUM  TURBULENT  VELOCITY  FLUCTUATIONS  AS  A  FUNCTION  OF  TIME 
AFTER  A  CONSTANT  COOLING  RATE  OF  1  ,8°c/H0UR  IS  APPLIED  TO  THE 
SURFACE. 


Constant  Cooling  Rate  :  Time  =  18  Hours 


Velocity  (M/Sec) 


FIGURE  4.6:  VELOCITY  DISTRIBUTIONS  OBTAINED  6  HOURS  AFTER  THE  APPLICATION  OF 
THE  CONSTANT  COOLING  RATE. 
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0  I - i—.rr— — I - 1 - 1 - 1 - 1 _ i _ *  ■ _ I _ I _ 

-0  -6  -4  -2  0  2  4 

Potential  Temperature 


FIGURE  4.7:  TEMPERATURE  DISTRIBUTION  OBTAINED  6  HOURS  AFTER  THE  APPLICATION 
OF  THE  CONSTANT  COOLING  RATE. 
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Friction  Velocity  (4 


Diurnal  Day  Variation 


Time  (  Hours  ) 

FIGURE  4.8:  DIURNAL  VARIATION  OF  THE  SURFACE  FRICTION  VELOCITY. 
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Surface  Temperature 


Diurnal  Day  Variation 


FIGURE  4.9: 


DIURNAL  VARIATION  OF  THE  SURFACE  TEMPERATURE. 


TurtiuJence  {  M/Sec  } 


Diurnal  Day  Variation 


FIGURE  4.10:  DIURNAL  VARIATION  OF  THE  MAXIMUM  TURBULENT  VELOCITY  VARIATION. 
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Diurnal  Day  Variation 


FIGURE  4.11:  DIURNAL  VARIATION  OF  THE  MAXIMUM  TURBULENT  LENGTH  SCALE. 
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Diurnal  Day  Variation  :  Time  =  40  Hours 


FIGURE  4.12:  WIND  VELOCITY  DISTRIBUTIONS  DURING  THE  MORNING  TRANSITION. 
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Diurnal  Day  Variation  :  Time  =  40  Hours 


FIGURE  4.13:  POTENTIAL  TEMPERATURE  DISTRIBUTION  DURING  THE  MORNING  TRANSITION. 
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Diurnal  Day  Variation  :  Time  -  44  Hours 


Velocity  (  M/Sec  ) 

FIGURE  4.14;  WIND  VELOCITY  DISTRIBUTIONS  DURING  AFTERNOON  UNSTABLE  CONDITIONS. 
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Diurnal  Day  Variation  :  Time  =  44  Hours 


FIGURE  4.15:  POTENTIAL  TEMPERATURE  DISTRIBUTION  DURING  AFTERNOON  UNSTABLE 
CONDITIONS. 
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Diurnal  Day  Variation  :  Time  =  52  Hours 


FIGURE  4.16:  WIND  VELOCITY  DISTRIBUTIONS  SHORTLY  AFTER  TRANSITION  TO  NOCTURNAL 
CONDITIONS. 
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Diurnal  Day  Variation  :  Time  =  52  Hours 


FIGURE  4.17:  POTENTIAL  TEMPERATURE  DISTRIBUTIONS  SHORTLY  AFTER  TRANSITION 
TO  NOCTURNAL  CONDITIONS. 


4-21 


5.  INCORPORATION  OF  CUMULUS  PARAMETERIZATION 
INTO  THE  TURBULENT  FLUX  RELATIONS 

A  large  part  of  turbulent  transport  in  the  atmosphere  involves  cloud 
dynamics.  Therefore,  before  any  subgrid  flux  parameterization  can  have  any 
general  applicability  for  a  mesoscale  model  it  must  represent  these  dynamics, 
at  least,  in  some  approximately  valid  manner.  The  cumulus  analysis  of 
Lewellen,  Sykes  and  Oliver  (1983)  provides  the  basis  for  extending  the 
parameterization  of  section  2  to  include  cumulus  transport.  As  long  as  the 
turbulent  kinetic  energy  equation,  Equation  2.1,  is  written  in  terms  of  the 
virtual  potential  temperature,  it  remains  ’unchanged  in  the  presence  of  change 
of  phase  of  the  water  in  the  atmosphere.  The  strong'  influence  of  the  change 
of  phase  comes  in  the  equations  for  6y  and  w’e^.  In  the  absence  of  any 
radiation  flux  divergence  it  is 


(5.1) 


which  is  the 


conserved  buoyancy  quantity  as  shown  in  Equation  5.12  of  L.S.O. 


De^  =  _  3W6’ 
Dt  ”  3z 


3W  61 

- -  + 

3z 


er  r’w 


(5.2) 


If  we  assume  that  the  vertical  flux  of  0V  needed  for  the  T.K.E.  equation  may 
be  given  by  a  quasi-equilibrium  version  of  Equation  5.9  of  L.S.O.  then  it  can 
be  written  as: 


w 1  6  1  v  = 


(5.  3) 
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with 


1/3 

1  +  8.6  R 

(5.4) 

R 

si' 

°i*r 

bC 

1 

(5.5) 

H  62 

Expressions  for 

p  cLnd  r'w1  may  ba  obtained  from  Equations  5*26  and  5«29 

r 

-  1/2  [j  +  Erf (q  / V^2 )J 

(5.6) 

w'  r ' 

-  (°r/°x)  (a  w'h'  "  b  w’6v  ) 

(5.7) 

with 

°r 

-  -^exp(-Q2/2) 

(5.8) 

Q 

*  (  K  *  hs)/°A 

(5.9) 

°x2 

-  a2  h'2  -  2ab  h'6’  ♦  b2  e;2 

.(5.10) 

a 

-  1  -  0.61  +  (  u  ^  +  0.61^  rj  f3y  hs 

(5.11) 

b 

-  £  1  +  (  u  1  ~1 )  rj  S-phg/Tp 

(5.12) 

V 

*v  2 

.  1  +  —  hsBf 

(5.13) 

By 

-  L/RyT 

(5.14) 
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The  expressions  for  wThf,  h’h1. 


6v- 


and  6 


,2 

v 


analogous 


to  Equation  5-3  are: 


v'h* 


-  *ash 


3h 

3z 


(5.15) 


(5.16) 


(5.17) 


(5.18) 


Thus 


o 


2 

X 


3h 

9z 


(5.19) 


This  set  of  equations  provides  a  complete  set  of  relations  for  two-phase 

turbulent  flow.  However,  the  character  of  the  turbulence  is  likely  to  be 

quite  different  under  conditionally  unstable  conditions  then  it  is  under 
either  fully  stable  or  unstable  conditions.  Equation  5.3  shows  a  given 
virtual  potential  temperature  gradient  may  be  stable  in  the  absence  of  cloud 
and  unstable  in  the  presence  of  the  cloud.  Further,  Equations  5.6,  5.9  and 
5.10  show  that  r  will  be  larger  in  the  presence  of  strong  turbulence.  Thus, 
there  is  a  natural  feedback  which  will  drive  the  flow  to  divide  between 

regions  of  high  turbulence  and  regions  of  little  or  no  turbulence.  It  is 
necessary  to  make  some  adjustments  in  our  transport  modeling  to  account  for 
this  intermi ttency .  We  define  intermit tency,  w  ,  as  the  volume  fraction  of 
the  sub-grid  ensemble  space  occupied  by  the  turbulent  regions.  We  argue  in 
LSO  that  the  principle  effect  this  has  on  our  transport  modeling  is  to 

increase  the  character ist ic  velocity  q  in  the  modeled  terms  by  uT^  so  that 
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they  are  based  on  the  average  turbulent  velocity  in  the  turbulent  region 
rather  than  the  average  over  the  layer  ensemble  space.  We  must  determine  how 
this  should  affect  Equations  5.2  to  5.18  and  how  this  new  variable  u>  should  be 

determined. 

Let  us  address  the  determination  of  oj  first.  Rather  than  attempt  to 
generate  an  equation  from  either  a  conditional  moment  approach  or  probability 
distribution  function  approach  (Kollman,  1984),  we  will  tentatively  rely  on  a 
simple  time  scale  relationship  to  determine  intermi ttency.  Under  stable 
conditions,  it  has  been  established  both  theoretically  and  observationally 
that  the  scale  of  3“D  turbulence  should  not  exceed  q/2N.  Theoret ically ,  this 
limit  is  based  on  the  energy  required  for  an  eddy  to  turn  over  in  the  presence 
of  stable  stratification.  Observationally,  it  agrees  with  the  limits 
necessarily  imposed  on  the  scale  to  obtain  agreement  with  the  Monin~Obukhov 
similarity  functions  in  the  stable  surface  layer  (Lewellen,  1981).  If  the 
turbulence  is  instead  permitted  to  be  intermittent  then  the  scale  can  exceed 
the  limit  of  q/2N.  Thus,  in  regions  where  the  scale  is  observed  to  exceed 
q/2N  we  can  expect  the  intermi ttency  to  be  proportional  to  q/NA,  i.e. 

w  -  Kq/NA  (5.20) 

Thus,  a  consistent  determination  of  A  will  through  Equation  5.20  also  provide 
a  determination  of  u.  For  purposes  of  preliminary  tests  we  will  stick  with 
essentially  the  same  A  algorithms  as  given  in  2.2  and  2.3  except  the  stability 
limit  on  A  is  applied  with  R  'based  on  the  moist  lapse  rate,  36v/3z  -  BT. 
Thus,  A  will  be  equal  to  A*  in  the  central  region  of  the  boundary  layer  and 
determined  by  CdA/dz]  *  0.65  at  both  the  top  and  the  bottom  of  the  boundary 
layer.  In  the  regions  where  q/2ft  would  have  yielded  a  smaller  A  the  flow  is 
intermittent  as  given  by  Equation  5.20.  The  biggest  weakness  of  this  approach 
to  i nt ermi ttency  is  that  it  builds  on  the  ad-hoc  nature  of  A.  However,  this 
is  not  as  damaging  as  it  might  first  appear  since  the  combination  wA  appears 
together  in  all  the  modeled  terms  except  the  diffusion  term  where  the  relevant 
combination  is  A/u>.  Since  there  is  theoretical  justification  for  taking  wA  as 
proportional  to  q/N  in  the  stable  regions  in  which  we  are  interested  most  of 
the  uncertainty  is  tied  to  the  diffusion  terms.  Of  course,  since  the  net 


effect  of  intermittency  is  to  increase  the  effect  of  the  diffusion  terms  it 
does  make  the  modeling  of  this  term  more  important. 


With  an  algorithm  for  determining  id  at  our  disposal  we  can  address  how 
the  quasi-equilibrium  expressions  should  be  affected  by  id.  By  following  the 
rule  that  q  in  the  modeled  terms  should  be  increased  by  id  ,  the 
quasi-equilibrium  expressions  may  be  written  as 


D(q2/2) 

Dt 


uf  w! 


3u 

3z 


V  f  w 1 


IX  +  £ 

3z  T 


we; 
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v_c± 

(D  3z 


(5.21) 
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R  =  min 


f!!v 
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0.25 
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-1/2 


(5.22) 


(5.23) 


(5. 2D) 


(5.25) 


(5.26) 


(5.27) 


(5.28) 
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It  is  also  desirable  to  adjust  the  cloud  functions  to  account  for  the 
nonguassian  distribution  of  the  humidity  and  temperature  fluctuations. 
Bougeault,  1981  has  shown  that  skewness  has  a  large  influence  on  the  r  and  or 
functions,  particularly  for  Q<  -  1 .  We  can  roughly  compensate  for  this  by 
allowing  the  guassian  expressions  given  in  Equation  5.6-5.10  to  apply  only  in 
the  turbulent  part  of  the  ensemble  volume.  The  expressions  then  become: 

1  ♦  Erf  (  Q /V2) 

°r  “  7f  exP  (  -Q2  /  2) 


(5.29) 


(5.30) 


Q..(h-hs)/ox  (5.31) 
with  given  by  Equation  5.19. 

If  we  wish  the  parameterization  to  include  precipitating  cumulus,  it  will 
also  be  necessary  to  add  a  term  to  the  w'h’  expression  which  represents  this. 
We  would  expect  this  to  be  accomplished  by  combining  a  sedimentation  rate 
which  is  a  function  of  droplet  size  Gunn  and  Kinzer  (1949.  with  a 
Marshall-Palmer  droplet  radius  distribution  as  a  function  of  liquid  water 
content.  The  liquid  water  content  can  be  obtained  following  the  same  analysis 
as  led  to  ? ,  or,  and  or  This  is  a  refinement  which  can  be  added  after  it  is 
demonstrated  that  the  basic  non-precipitating  parameterization  is  valid. 


As  a  test  of  our  algorithm  for  determining  intermittency ,  we  have 
recomputed  the  convective  boundary  layer  of  Section  4  using 
Equations  5.20-5.27.  Without  including  humidity  we  can  compare  with  the 
intermittency  observations  of  Deardorff,  Willis  and  Stockton  J98o, 
capping  stable  layer.  The  potential  temperature  profiles  of  Figure  5.1  and 
the  heat  flux  profiles  of  Figure  5.2  are  modified  very  little  from  those 
presented  earlier  for  the  Q.E.  result  of  Section  4.  The  intermittency 

results  in  Figure  5.3  are  in  quite  reasonable  agreement  with  Deardorff, 
et.al.’s  observation  points  for  K-1  in  Equation  5.20.  Although 
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Free  Convection 


FIGURE  5.1:  COMPARISON  OF  THE  TEMPERATURE  DISTRIBUTION  PREDICTED  FOR  FREE 

CONVECTION  BY  THE  INTERMITTENT  MODEL  WITH  THAT  FOR  THE  STANDARD 
Q.  E.  MODEL. 
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Heat  Flux 

FIGURE  5.2:  COMPARISONS  OF  THE  HEAT  FLUX  DISTRIBUTIONS  WHICH  GO  WITH  THE 
TEMPERATURE  DISTRIBUTIONS  OF  FIGURE  5.1. 
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Free  Convection 


FIGURE  5.3:  MODEL  PREDICTIONS  OF  THE  INTERMITTENCY  WITHIN  THE  INTERFACIAL  LAYER 
AT  THE  TOP  OF  A  FREE  CONVECTION  LAYER  AS  COMPARED  WITH  THE  DATA 
OF  DEARDORFF,  ET.AL  (1980).  ZETA  IS  THE  DISTANCE  FROM  THE  HEIGHT 
AT  WHICH  THE  HEAT  FLUX  GOES  TO  ZERO  NORMALIZED  BY  THE  INTERFACIAL 
THICKNESS. 
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intermittency  does  not  play  a  strong  dynamical  role  in  this  flow,  we  take  it 
as  a  very  encouraging  sign  that  our  intermittency  algorithm  is  compatibxe  with 
the  observations. 

In  checking  this  model  for  the  conditions  presented  in  LSO,  it  appears 
that  neglecting  terms  like  w'w'r'  in  the  w7 e;  equation  may  not  be  acceptable, 
although  this  would  seem  appropriate  for  a  quasi- equilibrium  approximation. 
This  type  term  can  be  included  in  a  straightforward  fashion  if  we  accept  the 
model  proposed  in  Equation  5.44  of  LSO.  In  this  case  it  is  only  necessary  to 
replace  r  in  Equations  5.24,  5.27  &  5.28  with 

r  +  aror/iD 

An  appropriate  value  of  a^,  will  need  to  be  determined  by  numerical 
experimentation.  A  value  of  unity  leads  to  qualitative  consistency  for  the 
example  considered  in  LSO. 
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6.  CONCLUDING  REMARKS 


We  have  presented  recommendations  for  a  subgrid  flux  parameterization 
scheme  appropriate  for  mesoscale  meteorological  models.  Initial  tests  of  this 
scheme,  as  compared  with  results  from  a  full  second-order  closure,  1 -D  model 
indicate  reasonable  results  when  model  grid  resolution  is  adequate  to 
determine  the  height  of  the  unstable  boundary  layer.  When  resolution  is  not 
adequate  for  this  it  will  be  desirable  to  include  a  prognostic  equation  for 
the  boundary  layer  thickness  or  incorporate  parameterization  for  an  elevated 
subgrid  shear  layer. 

Our  scheme  incorporates  cumulus  parameterization  as  an  integral  part  of 
the  quasi-equilibrium  formulation.  Our  tests  suggest  that  intermittency,  a 
fundamental  feature  of  cumulus  cloud  environments,  may  be  represented  in  the 
model  in  a  relatively  simple  fashion.  These  initial  results  indicate  that 
inclusion  of  the  effects  of  cumulus  dynamics  into  •  the  quasi-equilibrium 
approximation  of  turbulent  transport  should  be  feasible.  We  recommend  further 
testing. 
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ABSTRACT 

The  parameterization  of  subgrid  scale  turbulent  fluxes  of  momentum, 
energy,  and  mass  in  mesoscale  meteorological  models  is  discussed.  Theoretical 
analyses  and  available  observational  data  are  used  to  describe  our  current 
understanding  of  the  interdependence  between  turbulent  diffusion  and  the 
resolved  scale  distributions  of  wind,  temperature,  and  species.  Emphasis  is 
placed  on  the  dependence  of  subgrid  parameterization  on  the  simulation 
resolution  in  ensemble  space  as  well  as  in  physical  space  and  time.  Specific 
recommendations  are  presented  for  the  particular  problem  of  dispersion  from 
point  sources.  (*) 


‘Presented  at  the  AMS  workshop  on  Urban  Boundary  Layers,  October,  1983, 
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1  •  INTRODUCTION 


The  turbulent  fluxes  of  mass,  momentum  and  energy  are  at  least  partially 
composed  of  motions  that  are  too  small  to  ever  hope  to  be  resolved  in  urban  or 
mesoscale  models.  These  turbulent  fluxes  control  the  interaction  of  the 
atmosphere  with  the  surface  and  the  dispersal  of  anything  released  within  the 
atmosphere.  Such  mesoscale  phenomena  as  the  sea-breeze  circulation, 
mountain-valley  circulations,  and  the  moisture  build-up  in  the  boundary  layer 
necessary  to  drive  convective  clouds  are  all  dominated  by  turbulent 

interactions  with  the  surface.  As  long  as  the  turbulent  motions  in  these 
interactions  cannot  be  resolved  in  the  simulation,  the  parameterization  of 
subgrid  scale  fluxes  will  play  an  important  role  in  determining  the  accuracy 
of  the  simulation  of  any  of  these  phenomena.  This  is  perhaps  most  evident  in 
the  simulations  of  the  dispersal  of  passive  tracers  in  the  atmosphere  where 
almost  all  the  motions  responsible  for  dispersion  are  unable  to  be  resolved  in 
the  meteorological  model  of  the  region  of  interest. 

In  addition  to  the  question  of  what  scale  of  motion  can  be  resolved  by  a. 

feasible  grid  system,  there  is  the  question  of  how  much  of  the  motion  we  wish 

to  resolve.  Flow  in  the  atmospheric  boundary  layer  inherently  contains  a 
turbulent  stochastic  component.  Even  if  one  were  able  to  accurately 

numerically  simulate  all  of  the  scales  of  motion  in  time  and  space  for  one 
particular  realization,  this  would  not  provide  a  precise  prediction  of  the 
motion  in  time  and  space  for  any  other  particular  realization.  In  general, 
what  we  would  like  to  simulate  is  the  ensemble  mean  flow  distribution  in  time 
and  space.  We  would  also  like  to  gain  some  information  about  the  variability 
of  particular  realizations  frcxn  this  mean.  Numerical  simulations  either 
involve  ensemble  averaging  of  the  equations  or  averaging  of  the  simulations 
realized.  For  nonhomogeneous ,  nonstationary  problems  the  choice  of  scales  to 
average  over  before  the  numerical  simulation  is  an  important  part  of  the 
problem.  The  larger  the  scales  over  which  the  equations  are  averaged  the  more 
uncertainty  introduced  by  the  closure  approximations,  but  the  smaller  the 
scales  the  larger  the  computational  requirements  and  the  more  averaging  which 
must  be  done  after  the  simulation  to  provide  proper  interpretat ion  of  the 
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results.  The  implications  of  the  choice  of  the  ensemble  averaging  scale  on 
subgrid  parameterizat ions  are  discussed  in  the  next  section. 

The  rest  of  the  paper  reviews  some  subgrid  parameterizat ion  schemes  and 
speculates  on  how  they  may  be  improved.  Particular  attention  is  given  to 
dispersion  from  point  sources  because  this  is  a  problem  where  accurate 
parameterization  of  the  unresolved  turbulence  is  essential,  and  because  this 
is  an  area  where  we  have  carried  our  speculations  further. 

2.  THE  INFLUENCE  OF  THE  SCALE  OF  ENSEM3LE  AVERAGING  ON  SUBGRID  FLUX 

PARAMETERIZATION 

The  ensemble  of  flows  of  interest  are  all  the  possible  flows  which 
satisfy  the  prescribed  input  data.  In  ideal  problems  this  input  data  may  be 
sufficient  to  relatively  tightly  constrain  this  ensemble  of  flows,  but  in 
attempts  to  simulate  local  meteorological  flows  which  occur  at  a  specific  time 
and  place  the  input  data  is  unlikely  to  provide  tight  constraints.  Such 
simulations  must  be  able  to  deal  with  relatively  large  variances  from  the 
resulting  ensemble  mean  solution.  This  may  be  accomplished  either  by  ensemble 
averaging  over  the  equations  prior  to  performing  the  numerical  simulations  or 
by  performing  an  average  over  a  number  of  different  simulations  which  fall 
within  the  uncertainty  of  the  prescribed  input  conditions.  A  different  subgrid 
scale  flux  paramet er ization  needs  to  be  used  for  these  two  approaches.  In 
practice,  the  distinction  between  the  two  approaches  is  usually  made  in  terms 
of  the  scales  of  motion.  In  so-called  large  eddy  simulations,  some  ensemble 
averaging  is  done  at  small  scales,  while  higher-order  closure  approaches  which 
deal  with  ensemble  averaged  equations  over  much  larger  scales  may  still 
resolve  some  of  the  dominant  eddies  in  the  flow. 

The  standard  ensemble  average  velocity  may  be  defined  as 


(A-2.1) 
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where  {2  is  the  set  of  all  velocity  fields  which  satisfy  the  prescribed  input 

conditions  and  p(a)  is  the  probability  density  function  of  any  particular 

velocity  field  indicated  by  a.  We  wish  to  divide  up  this  domain  ensemble  into 

a  set  of  sub  ensembles  so  that  Each  is  the  set  of  all  velocity 

i 

fields  which  satisfy  the  prescribed  input  conditions  and  are  similar  at  scales 
of  motion  greater  than  some  scale  £e.  A  sub  ensemble  average 


<u(x,t)>i  -  f  u(x_, t )p( a) da  (A-2.2) 

^i 

is  then  equivalent  to  a  particular  large  eddy  simulation  of  the  velocity. 

In  order  to  make  this  concept  more  explicit,  consider  the  velocity  field 
u(x,t)  existing  in  a  space  domain  V0  during  a  time  interval  T0.  As  long  as  the 
Reynolds  number  based  on  a  characteristic  velocity  and  a  characteristic  length 
in  the  domain  is  large,  u  at  any  x  and  t  can  assume  a  wide  variety  of  values 
with  the  probability  of  any  individual  continuous  field  indicated  by  a 
specified  by  the  probability  density  function  p(a).  From  this  variety  of 
possibilities  let  us  choose  one  particular  continuous  field.  Having  chosen 

this  unique  field  we  average,  or  filter,  it  over  an  arbitrarily  chosen  scale, 

£  to-  eliminate  all  fine  scale  motions  of  a  scale  less  than  £0.  It  should  then 
be  possible  to  chose  other  neighboring  fields  that  when  similarly  filtered 
vary  at  most  by  an  e  from  that  chosen  first.  All  such  fields  are  included  in 
the  set.  Q-. .  As  long  as  the  Reynolds  number  based  on  £e  and  the  characteristic 
velocity  in  £21  is  large,  their  will  also  be  a  large  number  of  individual 
distributions  which  fall  within  fii .  The  process  can  be  repeated  to  form  a  Q2 
set,  etc.  The  number  of  possible  sets  is  also  large  as  long  as  £e  is 
significantly  less  than  the  dominant  shear  layer  thickness  in  the 

computational  domain. 

The  choice  of  £e  determines  the  break  between  resolved  and  unresolved 
motions.  Since  a  given  grid  cannot  resolve  a  motion  which  is  of  a 

characteristic  length  less  than  about  3  times  that  of  the  grid,  this  sets  a 
lower  bound  on  £g.  The  upper  bound  is  the  characteristic  length  of  the  domain 


of  simulation.  The  lower  bound  is  appropriate  for  Large  Eddy  Simulations  (LES) 
where  the  goal  is  to  resolve  as  much  of  the  motion  as  possible.  The  upper 
bound  is  appropriate  when  the  goal  is  to  resolve  only  the  mean  motion  in  the 
time  and  space  domain. 

3-  A  HIERARCHY  OF  SUBGRID  CLOSURE  ASSUMPTIONS 

The  ensemble  averaged  equations  are  precisely  defined  for  any  ensemble  as 
the  Reynolds  averaged  Navier-Stokes  equations 
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The  equations  are  well-defined  but  not  closed.  Some  assumptions  relating 
higher-order  correlations  must  be  made  to  close  the  system. 

The  assumptions  made  to  close  the  ensemble  averaged  equations  control 
what  parameterizations  a"e  appropriate  for  subgrid  turbulent  fluxes.  It  should 
perhaps  be  noted  that  when  Equation  3*2  is  used  for  LES  by  choosing  lQ  smaller 
than  the  scale  of  the  large  eddies,  there  is  no  inclusion  of  Leonard  stress 
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terms  (Leonard,  197*0  as  there  is  when  the  LES  equations  are  derived  by 
spatially  filtering  the  equations  of  motion.  The  u^  obtained  from  the  ensemble 
average  over  the  scale  of  is  thus  somewhat  different  than  the  u^  obtained 
by  spatially  filtering  over  a  scale  of  £e.  However,  this  difference  is 
generally  masked  by  approximations  in  the  Reynolds  stress  closure. 

The  simplest  closure  with  some  general  applicability  is  Smagor insky  ’  s 
(1963)  model  which  replaces  Equation  3*3  with  the  approximation 


(A-3.4) 


where 


(A-3.5) 


and  the  isotropic  part  of  the  Reynolds  stress  is  absorbed  within  the  mean 
pressure. 

The  length  scale  A  should  be  an  order  one  fraction  of  £e.  If  £e,  and 
consequently  A,  are  chosen  as  small  compared  to  the  large  eddies  expected  in 
the  field,  then  the  resulting  simulation  provides  a  single  realization  of  the 
large  eddies  in  the  flow.  In  this  case  the  desired  domain  ensemble  average  is 
obtained  by  averaging  the  results  over  a  suitable  number  of  these  sub-ensemble 
flows  or  large  eddy  simulations.  These  individual  realizations  may  be  obtained 
by  varying  the  input  conditions  within  the  constraints  of  the  specified 
accuracy  of  the  input  conditions.  On  the  other  hand,  if  £e  is  chosen  larger 
than  the  large  eddies  expected  in  the  field,  and  A  is  appropriately  chosen, 
then  the  simulation  will  approximate  that  of  the  ensemble  average  flow  over 
the  domain  without  any  further  averaging  of  the  simulation  results.  Of  course, 
at  the  level  of  closure  of  Equation  3.^  one  should  expect  the  first  LES 
procedure  to  provide  a  much  more  .accurate  simulation  in  return  for  the  much 
larger  computational  effort  required  to  accomplish  it. 
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If  A  is  chosen  to  be  larger  than  the  mesh  grid  but  still  sufficiently 
small  that  the  eddies  to  be  represented  by  (3-^)  are  well  within  the  inertial 
subrange  of  the  turbulent  spectra,  then  the  domain  ensemble  average  over  LES 
sub-ensembles  should  agree  with  the  domain  ensemble  results  for  smaller  A, 

although  the  individual  sub-ensemble  flows  may  be  quite  different  as  a 

function  of  x  and  t.  If  the  initial  input  for  the  two  simulations  are 
identical  then  the  two  simulations  should  initially  differ  only  in  their  small 
scale  eddy  motions  but  this  should  spread  to  differences  in  the  larger  eddies 
as  simulation  time  evolves.  The  higher  resolution  results  should  be  superior 
if  the  assumptions  on  which  (3*^)  a. re  based  are  not  completely  valid.  The  high 
resolution  results  will  also  provide  a  superior  simulation  when  the  resolution 
of  the  input  conditions  is  higher  than  the  low  resolution  LES.  If  A  is  chosen 
sufficiently  small  that  the  eddies  to  be  represented  by  (3*^)  are  well  within 
the  inertial  subrange  of  the  turbulent  spectra  then  it  is  possible  to  derive  a 
theoretical  value  for  c3.  However,  researchers  using  this  approach  have  found 
it  more  satisfactory  to  use  a  value  which  is  significantly  less  than  this 
theoretical  value  (cQ  *  0.065)  (Schumann,  1975)*  Equation  (3*^0  appears  to 

provide  reasonable  results  when  ujuj  ^  uiuj*  This  imposes  a  stringent 
requirement  on  grid  lengths  close  to  a  surface.  In  order  to  overcome  this  LES 
simulations  have  either  made  modifications  to  (3-^)  close  to  surface  or 

applied  boundary  conditions  which  incorporated  turbulent  surface  layer 
relationships. 

Close  to  a  surface  when  ujuj  becomes  of  the  same  order  as  u^Uj,  the 
assumptions  on  which  Equation  3*^  she  based  are  violated,  and  it  is  desirable 
to  transition  to  a  form  appropriate  for  ensemble  averaging  over  a  scale  which 
is  larger  than  the  grid  size  used  to  resolve  the  mean  flow  close  to  the  wall. 
The  common  mixing  length  model  which  is  valid  for  the  mean  Reynolds  stresses 
in  the  near-wall  region  is  similar  to  Equation  3A  but  with  A  replaced  by  a 
mixing  length  proportional  to  the  distance  from  the  surface.  The  transition 
formula  suggested  by  Schumann  (1975)  is 
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where  comes  from  the  mixing  length  formula  near  the  wall  and  is  not  allowed 
to  exceed  some  upper  bound  depending  upon  the  gri'd  mesh.  The  angular  brackets 
denote  a  transverse  spatial  average  which  due  to  the  symmetry  of  his  problem 
is  roughly  equivalent  to  an  ensemble  average.  Close  to  the  wall, 
S  ,-<S^  •>  <S^>  appreciably  less  than  one  on  average  and  the  first  term  on 

the  R.H.S.  of  (3.6)  is  much  less  than  the  last  term,  while  far  from  the  wall 
<S^j>  0  and  Equation  3*6  approaches  Equation  3*4. 

When  Equation  3.4  is  viewed  as  an  approximation  to  Equation  3*3  it 
becomes  apparent  that  the  non  local  dependence  involved  in  the  advection  and 
diffusions  terms  has  been  omitted,  the  influence  of  stratification  has  been 
ignored  and  it  has  been  assumed  that  the  unresolved  turbulence  is  isotropic. 
Under  these  conditions  the  dissipation  term  appearing  in  (3*3)  maY  te 
approximated  as 


e  -  0(W3«  6 i j / A  (A-3.7) 

and  the  equation  for  u{u[-q2  from  Equation  3-3  is  consistent  with  Equation 
3.4.  More  faithful  representations  of  Equation  3*3  should  perm  t  ie  to  be 
increased,  with  resultant  computational  savings,  without  significant 
sacrifices  in  accuracy.  Of  course,  the  computational  savings  are  only 
realized,  if  the  additional  computations  required  for  Equation  3-3  do  not 
exceed  the  savings  accrued  by  relaxed  space,  time,  and  ensemble  resolution 

requirements. 

Second-order  closure  models  for  Equation  3-3  appropriate  for  variable 
density  flows  have  been  investigated  by  a  number  of  researchers.  Recent 
reviews  have  been  given  by  Lewellen  (1981),  Zeman  (1981),  and  Mellon  and 
lamada  (1982).  A  model  of  Equation  3-3  which  has  proved  satisfactory  for  a 
number  of  applications  is 
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Compatible  equations  for  u>6',  6’  and  A  which  may  be  found  in  Lewellen  (1981) 


are  required  to  close  the  set  at  this  level  for  the  general  case.  Along  with 
different  forms  of  the  modeled  terms  in  Equation  3-8,  a  number  of 
approximations  intermediate  to  Equations  3*4  or  3*8  have  been  tried  in  the 
literature.  A  somewhat  more  general  formula  than  (3*4)  is  to  carry  an  equation 

similar  to  the  contraction  of  Equation  3-8  for  the  unresolved  turbulent 

2  - 

kinetic  energy  q  /2  and  then  to  model  Uj_Uj  as. 


(A-3.9) 


This  turbulent-kinetic-energy-transport  model  for  closure  has  been  quite 
popular  for  simulations  which  have  ensemble  averaged  over  all  scales  of 
motion,  and  used  by  some  LES  models,  e.g.  Schumann  (1975),  and 
Deardorff  (1980). 

The  next  level  of  approximation  is  to  allow  for  anisotropic  behavior  of 
vT.  This  can  be  obtained  by  combining  the  full  equation  for  from  Equation 

3.8  with  an  algebraic  approximation  to  Equation  3-8  for  the  relationship 
between  the  individual  components  of  ujul.  This  "algebraic  Reynolds  stress" 
model  carries  more  flavor  of  second-order  closure  than  using  Equation  3-9 
alone.  We  have  called  this  a  quasi-equilibrium  approximation  to  Equation  3-3 
in  our  use  of  it  to  compute  the  diurnal  variation  in  the  planetary  boundary 
layer  (Lewellen,  et.al.,  1974) .  Mellor  and  Yamada  (1974)  term  it  their  level 
2  1/2  approximation  to  second-order  closure. 

When  the  boundary  layer  assumption  is  valid,  i.e.,  when  gradients  in  only 
one  direction  are  important  in  the  production  terms  of  Equation  3-8,  the 
"algebraic  Reynolds  stress"  can  be  written  compactly  similar  to  the  form  given 
by  Yamada  and  Mellor  (1975)  as 


(A-3.10) 
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This  formulation  provides  for  the  important  influence  of  buoyancy  on  the 
eddy  viscosities  in  Equations  3-10  and  3-11 »  and,  has  been  used  by  Yamada 
(1979,  1983)  in  a  number  of  3~D  simulations  of  atmospheric  boundary  layers. 

The  horizontal  wind  variances  needed  for  dispersal  estimates  may  be 
obtained  the  same  way,  but  a  better  estimate  is  obtained  by  considering 
Equation  3*8  to  refer  only  to  the  fully  3“D  small  scale  turbulence.  Estimates 
of  the  ratio  of  the  horizontal  length  scale  of  the  turbulence  to  the  vertical 
length  scale  can  then  be  used  to  adjust  the  horizontal  wind  variance  so  that: 


u 
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(A-3.15) 
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It  is  important  to  note  that  this  1-D  version  of  quasi-equilibrium  is 
appropriate  for  boundary  layer  problems  with  £e>>6  where  the  simulation 
corresponds  to  the  domain  ensemble  average.  If  this  closure  were  used  for  a 
LES  with  ie<6  then  the  1-D  approximation  used  in  the  algebraic  reduction  would 
be  violated  by  the  inherent  3“D  character  of  the  large  eddies  in  the  boundary 
layer.  On  the  other  hand,  a  fully  3~D  algebraic  Reynolds  stress  formulation  is 
not  much  simpler  than  dealing  directly  with  Equation  3*8. 

The  only  computation  attempted  for  a  fully  3~D,  unsteady  Reynolds  stress 
model  has  been  that  of  Deardorff  (1972a)  of  flow  in  the  planetary  boundary 
layer.  He  used  £g  of  the  same  order  as  the  grid  mesh  to  perform  a  LES.  In  his 
numerical  simulations  which  imposed  a  logarithmic  law-of-the-wall  at  the 
surface  rather-  than  attempt  to  resolve  the  surface  layer,  the  added 
sophistication  of  using  a  modeled  Equation  3*3  did  not  provide  significantly 
different  results  from  that  obtained  using  the  simpler  Equation  3*1*.  We  would 
not  conclude,  as  some  have,  from  this  single  experiment  that  higher-order 
closure  is  not  useful  in  any  LES.  The  extent  to  which  Equation  3.3  can  be 
closed  to  faithfully  represent  more  of  the  turbulent  spectrum,  will  eventually 
determine  the  extent  to  which  second  order  modeling  can  make  LES  more 
practical. 

We  believe  the  marriage  of  LES  and  second-order  closure  (SOC)  should 
provide  benefits  unattainable  by  either  at  the  present  time.  The  relaxation  of 
the  grid  resolution  problems  at  the  higher  frequency  end  of  the  spectrum  of 
motions  should  allow  the  LES  to  concentrate  on  the  eddies  which  dominate 
turbulent  transfer  in  a  more  efficient  manner.  Conversely,  permitting  the  SOC 
model  to  directly  resolve  the  largest  eddies  should  provide  a  generality  to 
SOC  which  has  so  far  eluded  research  modelers.  The  2-D  LES  of  Lewellen,  Teske 
and  Sheng  (1980)  for  roll  vortices  in  the  planetary  boundary  layer  and  by 
Sykes  and  Lewellen  (1982)  of  Kelvin-Helmholtz  wave  breaking  in  a  shear  layer 
have  been  a  step  in  this  direction.  These  simulations  had  a  stronger  than 
desired  dependence  on  since  the  SOC  model  was  not  general  enough  to 

precisely  represent  the  effect  of  large  2-D  eddies  and  any  large  3-D.  eddies 
could  not  be  represented  by  the  2-D  simulation,  but  judicious  choices  of  £e 
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permitted  very  reasonable  simulations  of  the  two  phenomena. 

In  most  regional  scale  models  it  will  be  necessary  to  set  ie  as  of  the 
order  of  1  km  or  larger  since  any  finer  mesh  than  a  few  hundred  meters  in  the 
horizontal  appears  to  impose  prohibitively  large  computational  requirements. 
This  means  that  the  boundary-level  eddies  must  be  ensemble  averaged  over  and 
the  straightforward  application  of  (3.^)  is  only  appropriate  at  most  for 
parameterizing  the  horizontal  fluxes  of  momentum.  The  vertical  flux  of 
momentum  will  be  controlled  by  the  next  to  last  term  in  Equation  3.6.  In  this 
case  the  most  sensitive  parameter  in  the  formulation  is  the  bounds  placed  on 
in  the  outer  part  of  the  boundary  layer.  The  most  popular  formulation  for 
this  is  O'Brien's  (1970)  which  prescribes  a  cubic  polynomial  fit  to  vf  by 
matching  value  and  slope  to  surface  layer  relationships,  at  some  assumed  top 
of  the  surface  layer  and  fixing  the  value  at  the  top  of  the  boundary  layer. 
This  formulation  depends  explicitly  on  a  separate  dynamic  prediction  of  the 
boundary  layer  height.  A  number  of  other  formulations  for  vf  have  been 
reviewed  by  Blackadar  (1979). 

We  believe  a  reasonable  compromise  for  ie  of  order  10  km  is  to  carry  a 

2 

dynamic  equation  for  q  obtained  from  modeling  Equation  3*3.  A  proper  model  of 
this  equation  on  this  scale  should  recognize  the  disparate  nature  of  the 
largely  2-D  mesoscale  eddies  from  the  3-D  eddies  equal  to  or  smaller  than  the 
boundary  layer  scale.  A  1 -D  algebraic  Reynolds  stress  model  may  then  be  used 
for  the  vertical  fluxes  of  momentum,  energy,  and  mass  while  an  approximation 
similar  to  Equation  3*9  is  used  for  the  horizontal  fluxes.  The  determination 
of  the  vertical  distribution  of  the  length  scale  for  the  1-D  algebraic 
Reynolds  stress  model  will  require  an  estimate  of  the  boundary  layer 
thickness.  This  requirement  is  compatible  with  the  recommendation  that 
integral  equations  for  momentum,  thermal,  and  turbulent  kinetic  energy 
boundary  layer  thickness  be  used  in  connection  with  the  surface  boundary 
conditions  discussed  in  the  next  section. 
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4.  SURFACE  BOUNDARY  CONDITION  CONSIDERATIONS 


Strong  gradients  occur  in  most  of  the  meteorological  variables  in  the 
immediate  vicinity  of  the  surface*  The  surface  boundary  conditions  which  are 
most  appropriate  for  any  numerical  simulation  then  depend  critically  on  the 
vertical  grid  resolution  next  to  the  surface*  The  straightf orward  use  of 
no-slip  conditions  and  no  flow  thru  the  solid  boundary  are  generally  not 
possible  due  to  inadequate  resolution  next  to  the  surface*  It  is  possible  to 
show  that  one  should  have  Az<<zQ  before  no-slip  conditions  are  appropriate  at 
z-zQ*  Since  this  is  not  feasible,  one  must  either  asymptotically  match  to  the 
surface  layer  conditions  at  some  greater  height  where  Az<<z  is  feasible  or 
integrate  across  the  surface  layer  to  obtain  appropriate  boundary  conditions* 

A  reasonable  choice  is  to  define  the  values  of  turbulent  fluxes  at  the 
lowest  point  in  terms  of  a  transport  coefficient,  i.e. 


(A-U.1) 
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where  specification  of  c^  and  Cq  is  part  of  the  subgrid  par ameteri zation.  When 


the  lower  grid  point  falls  within  the  surface  layer  these  coefficients  may  be 


obtained  by  integrating  over  the  Monin-Obukhov  functions  (Paulson,  1970). 
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This  procedure  breaks  down  when  either  the  turbulence  within  the  surface 
layer  is  not  completely  controlled  locally  within  the  surface  layer,  or  the 
bottom  grid  point  falls  outside  the  surface  layer.  Deardorff  (1972a)  has 
provided  algorithms  to  circumvent  both  these  limitations.  The  first  is 
accomplished  by  replacing  u#  by  0.7w*  when  this  fraction  of  the  characteristic 
convective  velocity  exceeds  u*.  Since  w#  depends  on  the  boundary  layer 
thickness  this  provides  the  non-local  determination  of  the  turbulent- 
interaction.  The  second  limitation  is  bvercome  by  assuming  shape  functions  of 
5/L  for  the  velocity  and  temperature  profiles  over  the  outer  97 . 5%  of  the 
boundary  layer.  This  allows  him  to  generalize  the  surface  layer  c^  and  Cq 
relationships  to  be  a  function  of  6/zq,  the  ratio  of  the  boundary  layer 
thickness  to  the  surface  roughness,  and 


Ri 
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a  bulk  Richardson  number.  These  relationships  may  then  be  applied  even  when  z1 
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greatly  exceeds  the  surface  layer  thickness.  The  overall  boundary-layer 
thickness  is  estimated  from  an  integral  of  the  temperature  equation  with  upper 
bounds  applied  under  certain  neutral  and  stable  conditions.  Deardorff  intended 
his  PBL  parameterization  for  global  weather  models  but  it  has  also  been  used 
for  some  limited  area  models. 

We  believe  Deardorff’s  scheme  for  parameterizing  cf  and  c0  can  be 
improved  by  replacing  Equations  4.1  and  4.2  with 

u'w*  »  -Cf  q  u  (A-4.5) 

1  u 

v'w  »  ~£f  qv  ( A— 4 . 6 ) 


w’ 0 '  -  -ceq( Q^-Sg)  (A-4.7) 

where  q  is  obtained  from  the  dynamic  equation  for  the  turbulent  kinetic 
energy.  The  coefficients  defined  in  Equations  4.5  to  4.7  may  be  parameterized 
in  terms  of  4  different  boundary  layer  thicknesses.  Two  velocity  displacement 
thicknesses,  6U  and  6y,  the  thermal  thickness,  Se,  and  the  turbulent  kinetic 
energy  thickness  6q2,  may  be  determined  by  integrating  the  momentum, 
temperature,  and  turbulent  kinetic  energy  equations.  We  plan  to  investigate 
such  a  scheme  in  our  work  for  the  Naval  .Environmental  Prediction  Research 
Facility. 

We  expect  the  appearance  of  q  in  the  surface  coefficient  definition  to 
directly  provide  the  more  global  dependence  of  the  coefficients  required  under 
unstable  conditions  and  hope  that  the  use  of  the  4  integral  constraints  will 
permit  more  of  the  influence  of  unsteady  and  various  horizontal  pressure 
gradient  effects  to  be  included  in  the  parameterization  than  is  possibxe  when 
only  the  thermal  boundary  layer  thickness  is  used. 
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5.  DISPERSION  FROM  POINT  SOURCES 

A  subgrid  scale  problem  of  considerable  importance  is  that  of  dispersion 
from  a  point  source  of  a  passive  tracer.  When  the  plume  is  small  relative  to  a 
few  grid  widths  of  the  meteorological  simulation  its  dispersion  is  completely 
controlled  by  the  subgrid  turbulence.  Due  to  the  importance  of  turbulent 
parameterization  for  this  phenomenon  and  because  we  have  recently  made  what  we 
believe  to  be  significant  progress  in  this  area,  we  will  discuss  this  in  some 
detail.  A  subgrid  plume  from  a  point  source  may  be  modeled  as  a  gaussian  plume 
with  the  spread  of  the  plume  obtained  from  a  dynamic  equation  derived  by 
averaging  over  the  dynamic  equation  for  the  turbulent  species  flux  (Lewellen, 
1981).  With  the  mean  wind  and  turbulence  quantities  considered  homogeneous 
within  the  subgrid  volume  this  yields 


(A-5.1 ) 


2 

where  op  is  the  spread  of  the  plume  perpendicular  to  its  mean  trajectory. 


Equation  5.1  provides  the  basis  for  using  pertinent  information  about  the 
wind  fluctuations,  the  normal  components  of  the  Reynolds  stress,  and  the 
length  scale  of  the  turbulence,  to  directly  determine  the  spread  of  a  subgrid 
plume  rather  than  depend  on  empirical  functions  of  downwind  distance  for 
different  stability  classes.  It  is  important  to  note  that  Equation  5.1 
provides  an  estimate  of  the  ensemble  plume  spread.  For  relatively  large  values 
of  A/q,  the  time  scale  of  the  turbulence,  a  large  part  of  this  spread  will  be 
in  ensemble  space  which  manifests  itself  as  a  meandering  plume,  i.e.,  an 
uncertainty  in  position  of  a  much  narrower  plume  in  physical  space.  The 
magnitude  of  this  uncertainty  in  concentration  may  be  analyzed  by  dealing 
directly  with  the  concentration  variance  equation.  When  this  is  also  assumed 
to  have  a  guassian  distribution  the  integrated  variance  equation  can  be 
written  simply  as  (Lewellen  and  Sykes,  1983) 


(A-5.2) 
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where  tc  is  the  time  scale  of  the  dissipation  of  o'  and  <(  y  denotes 
integration  over  the  cross  section  of  the  plume.  The  Guassian  distribution 
has  been  shown  to  be  a  reasonable  approximation  for  homogeneous  turbulence  by 
Sykes,  Lewellen  and  Parker  (1983)  and  the  predictions  obtained  from  (5.2)  have 
been  demonstrated  to  be  compatible  with  the  data  on  concentration  variances 
from  a  point  source  in  a  wind  tunnel  boundary  layer  taken  by  Fackrell  and 
Robins  (1982).  Both  the  data  and  the  analyses  show  that  oc/c,  the  ratio  of 
the  standard  deviation  of  the  concentration  fluctuations  to  the  mean  value, 
can  often  exceed  1  even  on  the  centerline  of  the  plume  as  shown  in  Figure  A.1 
and  can  be  much  larger  in  the  edges  of  the  plume. 

The  calculation  of  the  concentration  variance  by  Equation  5.2,  or  similar 
equation,  provides  information  that  cannot  be  readily  obtained  from  releasing 
random  particles  into  even  a  completely  known  flow  unless  sufficient  particles 
are  tracked  to  provide  some  resolution  in  ensemble  space  as  well  as  physical 
space. 

6.  THE  USE  OF  A  VARIANCE  CALCULATION  AS  A  MEASURE  OF  PREDICTABILITY 

An  estimate  of  the  ensemble  variance  of  the  velocity  field  provides 
important  complementary  information  to  the  estimate  of  the  mean.  It  provides 
a  measure  of  the  predictability  of  the  velocity  and  is  essential  to  providing 
simulations  of  dispersion  of  a  passive  tracer  in  the  flow.  Uncertainty  in  the 
initial  conditions  of  a  particular  simulation  may  be  used  to  specify  initial 
values  for  the  velocity  variance  as  a  function  of  space.  In  a  SOC  simulation 
the  evolution  of  this  velocity  variance  field  may  be  tracked  as  the  mean  field 
is  tracked.  Subsequent  ratios  of  (ou/u)f  as  a  function  of  x  and  t  provide  a 
measure  of  how  well  the  velocity  field  is  specified  for  the  particular  input 
constraints.  Initially  ou  is  determined  completely  by  the  uncertainty  in  the 
initial  conditions.  However,  at  subsequent  times  both  the  internal  production 
terms  of  Equation  3-3  and  uncertainties  in  the  boundary  conditions  will 
contribute  to  au*  Even  if  it  were  possible  to  very  closely  define  the  initial 
conditions  and  boundary  conditions  so  that  cd=0  both  initially  and  on  the 
boundaries,  the  production  terms  in  Equation  3*3  would  still  impose  a 
variance  on  the  field  which  may  be  expected  to  grow  with  time  until  some  type 
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equilibrium  between  production  and  dissipation  is  reached.  Added  precision  in 
specifying  the  initial  conditions  will  only  reduce  the  variability  in  the 
field  at  times  earlier  than  this  equilibrium  time.  The  faithful  modeling  of 
Equation  3.3  should  thus  provide  a  useful  means  of  investigating  the  impact  of 
various  modeling  improvements  on  predictability  as  a  function  of  time. 

When  l  is  set  at  less  than  the  characteristic  scale  of  the  domain  so 
e 

that  individual  eddies  are  simulated,  the  effect  of  uncertainties  in  input 
conditions  must  be  determined  by  sensitivity  runs.  The  variability  introduced 
by  an  inverse  cascading  of  energy  from  high  frequency  motions  to  low 
frequency,  what  Leslie  and  Quarini  (1978)  refer  to  as  subgrid  backscatter, 
will  be  more  difficult  to  assess.  It  might  be  possible  to  assess  this  by 
investigating  the  sensitivity  of  the  LES  to  random  perturbations  in  the  LES 
field  which  are  consistent  with  the  variance  computed  for  the  sub  £e  motions. 

The  velocity  variance  determines  the  dispersion  of  a  passive  tracer  in 
the  flow.  It  is  generally  recognized  that  the  high  frequency  velocity  motions 
are  responsible  for  the  diffusion  of  a  tracer  as  a  puff  of  tracer  material  is 
transported  by  the  low  frequency  velocity  motions.  What  is  not  so  well 
recognized  is  that  intermediate  range  motions  may  either  diffuse  the  puff  or 
introduce  an  uncertainty  in  the  transport  of  the  puff.  All  of  the  variance 
serves  to  spread  the  puff  in  ensemble  space  but  only  the  motions  of  a  scale 
less  than  the  scale  of  the  puff  serve  to  diffuse  the  instantaneous  puff.  The 
rest  of  the  ensemble  spread  manifests  itself  as  a  meander  of  the  plume. 

The  contributions  of  these  meanders  to  the  diffusions  of  a  time-averaged 
puff  are  determined  by  the  relative  time  scale  of  the  meanders  and  the 
sampling  time,  A  mesoscale  gap  in  the  wind  energy  spectrum  (van  der  Hoven, 
1957)  is  often  evoked  as  justification  for  a  distinct  division  between  mean, 
transporting  winds,  and  fluctuating,  diffusing  winds.  This  gap  is  generally 
not  as  pronounced  as  theoreticians  would  desire.  The  wind  energy  spectrum  from 
the  100m  level  of  a  tower  at  the  Kincaid  power  plant  site  in  Illinois  is  given 
in  Figure  A. 2  for  2-3  week  periods  in  1980  as  taken  from  Lewellen,  Sykes,  and 
Parker  (1983).  Some  interaction  boundaries  have  been  sketched  on  the  spectrum 
to  indicate  what  part  of  the  spectrum  may  be  expected  to  contribute  to 
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FREQUENCY  CCYCLE3/H0UR) 

T'TCURK  A. 2:  Spectrum  of  the  horizontal  wind  f  luctuat  Ions  nt  the  Kincaid  site  (from 

Lewe  1  1  on  et.al.  I<JH3)  with  the  dashed  lines.  A,  H,  C,  I),  K ,  and  K  imllc.it  hip. 
the  boundaries  of  different  interact  Ions  between  the  plume  and  the  atmospheric 
motions.  These  Interactions  are  discussed  In  I  he  text. 


different  effects. 


Line  A  is  intended  to  represent  the  bound  of  motions  included  in  the  mean 
transport  of  the  plume.  The  motions  at  frequencies  less  than  this  bound  must 
be  adequately  resolved  so  that  the  transport  of  a  plume  element  may  be 
tracked.  There  must  be  a  consistency  between  the  time  for  which  the  plume  is 
tracked  and  the  spatial  resolution  of  the  wind  motion.  In  Figure  A.  2  the 
transition  between  transport  and  turbulence  has  been  arbitrarily  place  at 
ft-0.5  cycle/hour.  This  would  permit  a  4  m/s  mean  wind  to  transport  the  plume 
approximately  30  km  before  spatial  correlation  between  the  transporting  wind 
at  the  tower  and  at  the  plume  is  lost.  The  shape  of  line  A  is  symbolic  of  the 
fact  that  there  is  not  a  discontinuous  break  at  this  bound,  but  rather  than 
there  is  a  transition  range  of  frequencies  over  which  the  ability  to  include 
the  wind  energy  into  the  mean  transport  is  lost. 

Line  B  represents  the  lower  frequency  bound  of  the  energy  which  can  be 
included  as  a  part  of  the  turbulence.  B  should  complement  A  in  such  a  way  as 
to  assure  that  there  is  no  source  or  sink  of  wind  kinetic  energy  in  this 
transition  from  mean  transport  winds  to  turbulent  winds. 

At  the  high  frequency  end  of  the  spectrum,  line  E  represents  the  lower 
bound  of  motions  which  contribute  directly  to  the  diffusion  of  the 
instantaneous  plume.  Line  F  represents  the  upper  bound  on  motions  which 
contribute  to  the  total  ensemble  variance  of  the  concentration.  Higher 
frequencies  contribute  to  the  dissipation  of  the  variance  rather  than  its 
production.  The  larger  the  spatial  spread  of  the  instantaneous  plume,  the 
lower  the  frequency  which  can  contribute  either  to  the  diffusion  of  this 
instantaneous  plume  or  the  dissipation  of  the  concentration  variance.  Thus, 
lines  E  and  F  will  move  to  the  left  to  lower  frequencies  as  the  plume  spreads 
downstream.  The  turbulent  kinetic  energy  between  B  and  E  is  responsible  for 
the  meander  of  the  instantaneous  plume.  If  the  plume  is  tracked  sufficiently 
far  downwind  of  the  stack  then  E  may  move  to  the  left  of  B.  If  the  additional 
spatial  information  is  available  so  that  the  A-B  boundary  still  correctly 
represents  the  transition  from  transport  to  turbulence,  then  when  E  crosses  to 
the  left  of  B  it  means  that  part  of  the  resolvable  transport  motion  is  now 
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contributing  to  the  "diffusion"  of  the  instantaneous  plume. 


The  meander  of  the  plume  driven  by  the  energy  between  boundaries  B  and  E 
can  contribute  to  either  the  time-averaged  diffusion  of  the  plume  or  to  the 
uncertainty  in  the  position  of  the  time-averaged  plume.  The  location  of  these 
boundaries  C  and  D  on  Figure  A. 2  are  determined  by  the  sampling  time  period. 
The  position  of  C  and  D  sketched  is  arbitrarily  set  at  a  sampling  frequency  of 
twice  per  hour.  Line  C  represents  the  bound  on  energy  affecting  the 
concentration  level  of  the  time-averaged  sample.  Motions  represented  by  energy 
to  the  left  of  boundary  C  move  the  time-averaged  plume  around  as  a  coherent 
entity  rather  than  contributing  to  its  diffusion.  Boundary  D  represents  the 
boundary  between  the  motions  which  contribute  to  the  time-averaged  variance  of 
the  concentration  and  that  which  contribute  to  the  time-averaged  diffusion. 
Energy  to  the  right  of  D  contributes  only  to  the  time-averaged  diffusion  of 
the  plume.  We  expect  the  contribution  of  energy  in  frequencies  greater  than  fs 
to  the  time-averaged  variance  to  fall  off  as  (fs/f)  approximately;  this 
determines  the  shape  of  D.  The  shape  C  is  harder  to  set  but  is  determined  by 
the  enhanced  diffusion  resulting  from  the  interaction  of  the  small  scale  inner 
plume  turbulence  with  the  distortions  of  the  plume  forced  by  the  large  scale 
motion.  As  sampling  time  is  reduced  lines  C  and  D  approach  lines  E  and  F, 
respectively. 

Our  purpose  here  is  not  to  try  to  precisely  define  the  shape  of  all  the 
boundaries  on  Figure  A. 2,  but  to  argue  that  such  boundaries  exist  and 
qualitatively  note  the  type  of  influence  the  energy  bounded  by  the  different 
lines  has  on  the  plume.  This  breakdown  of  the  interactions  of  different  scales 
of  motion  illustrates  that  a  level  of  uncertainty  is  an  inherent  part  of  plume 
dispersion  which  can  not  be  eliminated  even  by  a  perfect  model.  However, 
improved  models  should  be  able  to  provide  an  estimate  of  the  variance  along 
with  their  predictions  of  the  mean  concentration  distribution. 
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7.  RELATED  SUBGRID  PARAMETERIZATION  PROBLEMS 
1 )  Subgrid  Scale  Sources 

Not  only  is  it  necessary  to  model  the  subgrid  transport  of  mass, 
momentum,  and  energy,  but  it  is  also  often  necessary  to  model  subgrid  sources 
(or  sinks)  of  one  or  more  of  these  quantities.  Relatively  less  work  has  been 
de.voted  to  this  subject  than  to  the  modeling  of  the  turbulent  transport  terms 
although  in  urban  models  the  influence  is  likely  to  be  nearly  as  important. 
Subgrid  surface  sources  generally  have  been  incorporated  into  the  effective 
roughness  of  the  surface.  This  is  appropriate  as  long  as:  (1)  the  height  of 
the  source  is  well  below  the  first  vertical  grid  point  above  the  surface,  and 
(2)  the  sources  are  relatively  uniformly  distributed  between  the  horizontal 
grid  points.  When  this  is  not  the  case  a  more  accurate  formulation  should 
improve  the  simulation,  particularly  if  the  detailed  influence  of  these 
sources  are  of  interest.  Pielke  (1981)  has  proposed  that  subgrid  scale 
terrain  forcing  can  be  neglected  in  mesoscale  models  if  the  subgrid  scale 
terrain  height  variance  is  small  compared  with  the  grid  resolvable  terrain 
variance  and  has  investigated  the  Spectra  of  terrain  height  variance  over 
parts  of  Virginia  (Pielke  and  Kennedy,  1980)  and  Colorado  (Young  and  Pielke, 
1983)  to  see  what  limits  on  mesoscale  resolution  this  imposes  as  long  as 
subgrid  terrain  is  neglected.  The  conclusion  is  that  this  limits  the  maximum 
allowable  horizontal  grid  spacing  to  100m  in  Colorado  and  1  km  in  Virginia. 
Subgrid  terrain  forcing  is  probably  even  more  important  than  that  implied  by 
the  ration  of  the  integrated  subgrid  terrain  height'  variance  to .  the  integral 
of  the  resolved  terrain  variance  due  to  the  role  of  subgrid  terrain  changes  in 
separating  the  boundary  layer  flow.  However,  even  this  conservative  estimate 
indicates  that  such  forcing  can  not  be  neglected  in  hydrostatic  mesoscale 
simulations  which  have  a  lower  bound  on  the  horizontal  grid  spacing  of 
approximately  5  km.  This  problem  clearly  warrants  more  attention  than  it  has 
received. 
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2)  Unresolved  Shear  Layers  Away  from  the  Surface 


Another  subgrid  phenomena  which  merits  work  in  mesoscale  models  is  that 
of  dealing  with  local  regions  of  high  gradients  which  cannot  be  resolved  by 
the  grid  mesh.  It  might  be  argued  that  the  grid  should  be  set  so  that  all 
shear  layers  of  importance  to  the  phenomena  under  investigation  could  be 
resolved,  but  this  is  not  likely  to  always  be  feasible.  Thus,  subgrid  flux 
parameterization  may  be  improved  if  it  can  be  formulated  in  such  a  way  as  to 
permit  thin  shear  layers  to  exist  within  the  subgrid  space  without  being  wiped 
out.  The  top  of  the  surface  mixed  layer  or  the  remnants  of  what  was  once  the 
top  of  the  mixed  layer  often  falls  into  the  category  of  a  shear  layer  which 
needs  such  parameterization.  A  potential  technique  for  this  might  be  to 
compute  the  subgrid  Richardson  number  in  the  vertical  direction  and  limit  it 
to  not  exceed  a  critical  value  of  approximately  1/1),  i.e.,  require  that 


Ri 


ATAz 

[  M2  +  q2] 


<  i/ii 


(A-7.1) 


The  rationale  for  imposing  this  limit  is  the  recognition  that  if  Ri  is 
significantly  bigger  than  this,  turbulence  cannot  be  supported  within  the 
layer.  Conversely,  if  Ri  is  significantly  less  than  this  value  the  layer  is 
liKeiy  to  be  unstable  and  grow  in  thickness.  Thus  a  simple,  rough 
parameterization  of  a  turbulent  shear  layer  may  be  accomplished  by  setting 
Ri»1/4.  If  the  computed  value  exceeds  1/4  then  a  subgrid  shear  layer  with  a 
Ri-1/4  may  be  placed  within  the  grid.  The  shear  layer  would  be  joined  by 
slopes  of  AT_/Az_ ,  Au_/Az_  on  the  lower  side  and  AT+/Az+,  Aud+/Az+,  on  the 
upper  side.  The  subgrid  location  of  the  shear  layer  could  only  be  determined 
if  other  considerations  are  introduced.  Compatibility  with  the  computed 
boundary  layer  thickness  could  provide  such  a  constraint  in  the  case  of  the 
shear  layer  at  the  top  of  the  boundary  layer. 


3)  Cumulus  Paramet erization 


Cloud  dynamics  play  a  strong  part  in  the  interaction  of  the  surface  with 
the  atmosphere.  Cumulus  clouds  present  a  particular  problem  to  mesoscale 
simulations  because  £e  is  sufficiently  large  that  individual  clouds  can  not  be 
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resolved  and  the  dynamical  role  of  clouds  must  be  included  in  the  subgrid 
parameterization.  When  the  clouds  are  fairly  uniformly  distributed  over  the  £e 
domain  then  the  principal  import  of  the  clouds  is  the  influence  of  the  latent 
heat  release  on  the  temperature  equation  and  the  subsequent  effect  of  this  on 
the  buoyant  generation  of  vertical  motions.  An  added  complication  is 
introduced  when  the  cloud  instability  begins  to  control  the  turbulent  motion. 
This  natural  instability  of  the  cloud  motion  tends  to  concentrate  the  vertical 
motions  and  causes  the  turbulence  to  be  distributed  very  intermittently  in 
time  and  space.  The  usual  turbulence  closure  approximations  assume  some 
uniformity  of  the  turbulence  within  the  ensemble  domain.  Strong  intermittency 
can  seriously  invalidate  these  approximations.  Thus  a  serious  attempt  to 
incorporate  cumulus  dynamics  into  a  general  representation  of  subgrid 
turbulent  fluxes,  rather  than  rely  on  ad  hoc  parameterization  of  the  influence 
of  cumulus  convection  such  as  those  given  by  Kuo  (1974)  and  Arakawa  and 
Schubert  (1974),  must  find  a  reasonable  closure  approximation  for  intermittent 
turbulence.  We  have  speculated  (Lewellen,  et.al.,  1983)  that  this  may  crudely 
be  accomplished  by  allowing  the  turbulent  time  scale  to  be  proportional  to  an 
intermittency  variable.  This  permits  a  more  rapid  redistribution  and 
dissipation  of  turbulent  energy.  It  remains  to  be  seen  how  well  we  can 

accomplish  the  critical  step  of  determining  an  appropriate  equation  for  the 
intermittency  variable. 

8.  CONCLUDING  REMARKS 

We  wish  to  emphasize  the  importance  of  correctly  interpreting  variability 
within  ensemble  space  when  attempting  to  simulate  mesoscale  temporal  and 
spatial  distributions  of  meteorological  variables.  The  appropriate  choices  for 
subgrid  flux  parameterization  depend  strongly  on  the  domain  in  ensemble  space 
a  modeler  wishes  the  simulation  to  represent.  The  subgrid  flux 
parameterization  is  easier  when  the  ensemble  domain  is  tightly  constrained, 
but  more  averaging  of  particular  realizations  is  required  to  provide  a  correct 
interpretation  of  the  simulation  results.  We  have  attempted  to  provide  some 
guideline  for  flux  parameterization  in  general  and  have  recommended  a  specific 
framework  which  we  will  be  testing  for  the  Naval  Environmental  Prediction 
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Research  Facility.  We  have  also  made  specific  recommendations  for  the 
particular  problem  of  dispersion  from  point  sources. 


A-26 


References 


Arakawa,  A.,  and  W.H.  Schubert,  1974:  "Interaction  of  a  Cumulus  cloud  ensemble 
with  the  large-scale  environment,  Part  I,"  J.  Atmos.  Sci.,  31.  pp. 
674-701. 

Blackadar,  A.K.,  1979:  "High-resolution  models  of  the  planetary  .boundary 

layers,"  Advances  in  Envir.  Sci.  and  Eng.,  1_,  pp.  50-85. 

Deardorff,  J.W.,  1972a:  "Numerical  investigation  of  neutral  and  unstable 

planetary  boundary  layers,"  J.  Atmos.  Sci . ,  29 ,  pp.  91-115. 

Deardorff,  J.W.,  1972b:  "Parameterization  of  the  planetary  boundary  layer  for 
use  in  general  circulation  models,"  Monthly  Weather  Review,  100,  pp. 
93-106. 

Fackrell,  J.E.,  and  A.G.  Robins,  1982a:  "Concentration  fluctuation  and  fluxes 
in  plumes  from  point  sources  in  a  turbulent  boundary  layer,"  J.  Fluid 
Mech. ,  117,  pp.  1-26. 

Kuo,  H.L.,  1974:  "Further  studies  of  the  parameterization  of  the  influence  of 
Cumulus  convection  on  large-scale  flow,"  J.  Atmos.  Sci.,  31,  pp. 

.  1232-1240. 

Leonard,  A.,  1974:  "On  the  energy  cascade  in  large-eddy  simulations  of 
turbulent  fluid  flows,"  Advances  in  Geophysics  1 8 A ,  pp.  237-248. 

Leslie,  D.C.,  and  G.L.  Quarini,  1979:  "The  application  of  turbulence  theory  to 
the  formulation  of  subgrid  modeling  procedures,"  J.  of  Fluid  Mech.,  91, 
pp.  65-91 . 

Lewellen,  W.S.,  1981 :  "Modeling  the  lowest  1  Km  of  the  atmosphere,", 
AGARD-AG-267.  Available  from  NTIS. 

Lewellen,  W.S.,  and  R.I.  Sykes  as  a  measure  of  natural  uncertainty  in  observed 
concentration  samples,"  Proc.  AMS  6th  Symp.  Turbulence  and  Diffusion, 
Boston,  MA.  "  .  ’  ~~  ~ 

Lewellen,  W.S.,  R.I.  Sykes,  and  D.A.  Oliver,  1983:  "Further  developments  of 
the  A.R.A.P.  model  for  the  atmospheric  marine  environment,"  A.R.A.P. 
Report  No.  488. 

Lewellen,  W.S.,  R.I.  Sykes,  and  S.F.  Parker,  1983:  "Analysis  of  plume 
variability  based  on  select  periods  of  the  Kincaid  data  set,"  ARAP  Report 
495,  May  1983. 


A-27 


Lewellen,  W.S.,  M.E.  Teske,  and  C.  duP.  Donaldson,  1  974:  "Turbulence  model  of 
diurnal  variations  in  the  planetary  boundary  layer,"  Proc.  of  the  1974 
Heat  Trans,  and  Fluid  Mech.  Inst.,  (edited  by  L.R.  Davies  and  R.E. 
Wilson),  Stanford  Univ.  Press,  Stanford,  CA,  pp.  301-319. 

Lewellen,  W.S.,  M.E.  Teske,  and  Y.P.  Sheng,  1980:  "Micrometeorological 
applications  of  a  second-order  closure  model  of  turbulent  transport," 
Turbulent  Shear  Flows  2,  Selected  papers  from  the  Second  Int'l.  Symp.  on 
Turbulent  Shear  Flows,  Imperial  College  London,  July  2-4,  1979, 
Springer-Verlag,  (New  York),  pp.  366-378. 

Mellor,  G.L.,  and  T.  Yamada,  1974:.  "A  hierarchy  of  turbulence  closure  models 
for  planetary  boundary  layers,"  J .  Atmos .  Sci . ,  31 ,  PP-  1791-1806. 

Mellor,  G.L.,  and  T.  Yamada,  1982  "Development  of  a  turbulence  closure  model 
for  geophysical  fluid  problems,"  Reviews  of  Geophysics  and  Space  Physics, 
20,  pp.  851-875. 

O'Brien,  J.J.,  1970:  "A  note  of  the  vertical  structure  of  eddy  exchange 

coefficients  in  the  planetary  boundary  layer,"  J.  Atmos.  Sci.,  27,  Dp. 
1213-1215. 

Paulson,  C  A.  1970:  "The  mathematical  representation  of  wind  speed  and 
temperature  profiles  in  the  unstable  atmospheric  surface  layer,"  J.  Appl. 
Meteorology,  9,  pp.  857-861  . 

Pielke,  R.A.,  1981:  "Mesoscale  numerical  modeling,"  Adv.  Geophys. ,  23,  pp. 

1 85-344. 

Pielke,  R.A.,  and  E.  Kennedy,  1980-  "Mesoscale  terrain  features,"  Report  No. 
UVA-ENV  SCI-MES0-1 980-1 ,  University  of  Virginia.  (Avail,  from  R.A. 
Pielke,  Dept,  of  Atmospheric  Science,  Colorado  State  Univ.,  Fort  Collins, 
Colorado  80523.) 

Schumann,  U.,  1975:  "Subgrid  scale  model  for  finite  difference  simulations  of 
turbulent  flow  in  plane  channels  and  annuli,"  J.  Comp.  Phys.,  18,  pp. 
376-404. 

Smagorinsky,  J.S.,  1963=  "General  circulation  experiments  with  the  primitive 
equations,  I:  The  basic  experiment,"  Monthly  Weather  Review,  91 ,  pp.  99- 

Sykes,  R.I.,  and  W.S.  Lewellen,  1982:  "A  numerical  study  of  breaking 
Kelvin-Helmholtz  billows  using  a  Reynolds-stress  turbulence  closure 
model,"  J.  Atmos.  Sci . ,  39,  pp.  1506-1520. 

Sykes,  R.I,  W.S.  Lewellen,  and  S.F.  Parker,  1984:  "A  turbulent  transport  model 
for  concentration  fluctuations  and  fluxes,"  J.  Fluid  Mech.,  139, 
pp.  193-218. 


A-28 


Van  der  Hoven,  I.,  1957:  "Power  spectrum  of  horizontal  wind  speed  in  the 
frequency  range  from  0.0007  to  900  cycles  per  hour,"  J .  Meteor . ,  1 4 ,  pp. 
1 60-169. 

Yamada,  T. ,  1979:  "An  application  of  a  three-dimensional  simplified 

second-moment  closure  model  to  study  atmospheric  effects  of  a  large 
cooling-pond,"  J.  Atmospheric  Environment,  pp.  693-704. 

Yamada,  T. ,  1983:  "Numerical  simulation  of  winds  and  turbulence  in  the 
California  geysers  area,"  Proc.  AMS,  6th  Symp.  Turbulence  and  Diffusion., 

Yamada,  T. ,  and  G.L.  Mellor,  1975:  "A  simulation  of  the  wangara  atmospheric 
boundary  layer  data,"  J.  Atmos.  Sci.,  32 ,  pp.  2309-2329. 

Young,  G.S.,  and  R.A.  Pielke,  1983:  "Application  of  terrain  height  variance 
spectra  to  mesoscale  modeling,"  submitted  to  JAS. 

Zeman,  0.,  1981:  "Progress  in  the  modeling  of  planetary  boundary  layers," 

Annual  Review  of  Fluid  Mech.,  13,  PP-  253~272. 


A-29 


APPENDIX  B 


PARAMETERIZATION  OF  THE  SURFACE  LAYER  UNDER  CONDITIONS 
APPROACHING  FREE  CONVECTION 


Monin-Obukhov  surface  similarity  contains  a  singularity  under  free 

convection  conditions.  We  wish  to  explore  the  parameterization  of  the  surface 

layer  under  conditions  approaching  this  singularity.  It  is  often  assumed 

1  /3 

under  these  conditions  that  ow  -  w#z  and  -  w*.  Indeed,  available  data 
taken  in  the  surface  layer  appears  to  support  this  limit  reasonably  well. 
However,  this  data  is  generally  taken  at  moderately  large  values  of  z  to  be 
sure  that  (-z/L)  is  very  large.  This  simple  form  cannot  hold  at  sufficiently 
small  z  because  must  be  reduced  to  zero  in  some  neighborhood  of  zQ.  Thus, 
we  expect  some  inner  region  of  the  free  convection  surface  layer  to  be 
dominated  by  a  cascading  process  induced  by  the  unsteady  shear  stress  produced 
by  the  horizontal  fluctuations  interacting  with  the  surface  roughness 
elements.  We  will  provide  an  estimate  of  the  height  of  this  inner  layer  and 
discuss  the  implications  of  this  layer  for  the  parameterization  of  surface 
boundary  conditions. 
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ANALYSIS 


In  order  to  provide  a  quantitative  estimate  of  this  inner  layer 
thickness,  let  us  approximate  the  time-averaged,  small-scale,  turbulent 
kinetic  energy  equation,  vertically  integrated  over  this  layer  as: 


0 


(B.1) 


The  first  term  represents  the  shear  production  of  high  frequency  energy 
resulting  from  the  relatively  low  frequency  horizontal  velocity  fluctuations 
interacting  with  the  surface.  The  second  term  represents  the  buoyant 
production  term  due  to  the  heat  flux  through  this  layer.  The  final  term 
represents  the  dissipation  of  high  frequency  energy  in  this  region  were  the 
dissipation  length  scale  is  assumed  to  be  inversely  proportional  to  the 
distance  from  the  surface.  The  coefficient  0.2  comes  from  the  ratio  of  the 
two  coefficients  b/a  in  the  standard  A.R.A.P.  model. 

We  have  defined  the  inner  layer  as  that  region  within  which  the  shear 
production  exceeds  the  buoyant  production.  Let  us  arbitrarily  define  6  as 
precisely  the  height  where  the  first  2  terms  of  Equation  1  become  equal. 
Above  this  height,  the  first  term  will  remain  essentially  constant  while  the 
second  term  continues  to  grow  linearly.  From  both  laboratory  (Willis  & 

Deardorff,  1974)  and  field  observations  (Panofsky,  et.al.,  1977),  the  u^ 
scaling  is  given  by 


(B.  2) 
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thus 
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Equations  (3)  and  (5)  give 
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Under  neutral  surface  layer  conditions  as  defined  by 


would  be  given  by 
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cf  -  0.16/ln  (6i/z0)  (B. 8) 

When  Equation  8  is  increased  by  50?  to  allow  for  some  instability  in  this 
layer,  and  used  for  c^.  in  Equation  6,  the  final  result  is 


Si  (ln  5i/zo)2  -  °->7  » 


(B.9) 


When  h  is  approximately  1km  and  zQ  is  approximately  0.01m,  Equation  9 
will  give  6i  »  4.6m.  Over  water  where  the  unstable  mixed  layer  may  be  as 

-4 

small  as  100m  and  zQ  as  small  as  1 0  m,  the  inner  surface  layer  thickness  will 
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be  as  small  as  =  0.3m. 


When  zQ  is  greater  than  approximately  0.1,  a  significant  part  of  this 
inner  layer  is  likely  to  be  occupied  by  the  canopy  layer  which  extends  out 
past  10zc  and  the  horizontal  velocity  fluctuations  will  play  an  even  more 
active  role  in  the  interaction  between  the  surface  and  the  atmosphere. 

Deardorff  .(1972)  in  his  parameterization  of  the  planetary  boundary  layer 
recognized  that  free  convection  conditions  can  not  extend  all  the  way  to  the 
surface.  His  reasonable  fix  for  this  condition  was  to  not  allow  cu  to  exceed 


2  c,,  or  cQ  to  exceed  3-33  cQ 

uN  e  eN 


where  he  defined  cu  and  Cq  as: 


(B. 1 0) 


c. 


u 


-  w'e* 

0 


CB. 11) 


and  c.,  and  cQ  are  the  neutral  values  of  these  transport  coefficients. 

UN  °N 

The  accuracy  of  cu  is  not  very  critical  in  this  limit  where  the  mean  wind 
is  approaching  zero  but  the  heat  flux  cannot  be  allowed  to  go  to  zero  also. 
Deardorff  imposes  a  lower  bound  of 


w'0*  free  convection  »  0.0019  (  0S  -  0m  'j  m.dg/sec  ( B .  1  2 )  • 


Our  inner  layer  analysis  presented  here  suggest  that  the  lower  bound  on 


is  given  by 


(B. 1 3) 


so  that 


we'f 


.  c. 


10  6i  y/2 

m  6i/zoy/ 


(B. 15) 


Consistent  with  Equation  9,  we  estimate  that  c#>  should  be  approximately  1.5 

1  9 

times  that  appropriate  for  a  neutral  surface  layer  of  thickness  6^,  i.e., 

cf^  =  0.33 /  (in  6i/z0j  (B.16) 

thus 

1/2(es  -  em  )  3/2  (B*17) 

For  typical  values  of  g/TQ  =  1/30.  =  Ik®.  an^  z0  =  0.01m, 

this  gives 

-  °-006  (»s  -  6m)  3/2 

which  is  approximately  three  times  that  given  by  Equation  12  when 
^or’  zo  =  0.1m  and  0S  -  6m  =  10  the  ratio  would  be  approximately 
10.  Thus,  for  large  temperature  differences  and  rough  surfaces  Equation  12 
can  underestimate  the  free  convection  heat  flux  by  as  much  as  an  order  of 
magnitude. 

Louis  (1979,  1982)  provides  a  PBL  parameterization  which  approaches  the 
free  convection  limit  with  heat  flux  given  by 


w'0'f 


.  c. 


(B.19) 


This  is  the  same  temperature  dependence  as  given  in  Equation  17  and  although 
Equation  19  is  independent  of  h  and  a  quite  different  dependence  on  z0,  it 
leads  to  an  expression  approximately  40?  less'  than  Equation  18  for  the  same 
conditions  as  specified  for  that  equation.  If  zQ  is  decreased  from  10  m  to 
10  ^m  the  coefficient  in  Equation  18  would  be  reduced  by  approximately  a 
factor  of  4  while  Equation  19  would  lead  to  a  factor  of  10  reduction. 


B-5 


CONCLUSIONS 


The  free  convection  surface  layer  conditions  of  ow  “  z  ^  will  not  extend 
to  the  surface  at  z  .  When  the  Monin-Obukhov  is  negative  and  close  to  zero 
the  convective  inner  layer  described  here  will  extend  from  z  *  20  z0  out  to 
z  •  where  is  given  by  Equation  9.  This  layer  will  not  occur 
when  -L  >  6„  *  . 

C  •  i.  « 
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APPENDIX  C 


TWO-DIMENSIONAL  TURBULENT  SIMULATIONS 

The  justification  for  two-dimensional  studies  comes  from  atmospheric 
observations,  and  this  allows  us  to  make  calculations  of  the  large  eddies 
explicitly  and  obtain  results  which  are  less  sensitive  to  closure  assumptions. 
We  have  utilized  two  different  models  for  these  studies;  one  is  the  A.R.A.P. 
closure  model  as  reported  by  Teske  and  Lewellen  (1979),  and  the  second  is  the 
mixing-length  parameterization  model  of  Mason  and  Sykes  (1982).  Both  models 
were  previously  used  to  study  boundary-layer  turbulence- profiles  and  employed 
artifices  which  reduced  the  reliability  of  the  entrainment  predictions.  Teske 
and  Lewellen  prescribed  numerical  damping  at  a  short  distance  above  the 

inversion,  while  Mason  and  Sykes  included  an  artificial  cooling  term  to  obtain 
a  steady  mean  state  at  the  expense  of  a  specific  relation  between  the  heat 
flux  and  mean  temperature  profiles.  In  the  studies  reported  here,  we  consider 
a  growing  boundary  layer  beneath  a  deep,  stably-stratified  layer  so  that 
internal  gravity  waves  are  generated  and  transport  momentum  and  energy  in  the 
vertical.  An  absorbing  layer  is  placed  at  the  upper  levels  to  prevent 

reflection  there. 

The  need  to  consider  a  growing  boundary  layer  forces  a  compromise  on  the 

length  of  the  numerical  integration  time;  we  wish  to  integrate  for  a  long 

time  so  that  initial  conditions  are  not  influencing  the  results,  but  we  do  not 
want  the  boundary  layer  to  grow  so  deep  as  to  make  the  computational  mesh 
inappropriate.  These  requirements  led  us  to  choose  the  following  parameters 
for  our  runs:  geostrophic  wind  speed,  Ug  -  10ms”1;  initial  boundary  layer 
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depth,  z^  -  500m;  initial  surface  heat  flow,  HQ  =  0.1  °Kms  1  ;  overlying 
stable  temperature  gradient,  r  -  6x10-^  °Km_1 .  The  latter  is  a  relatively 
strong  gradient,  but  these  parameters  produce  a  boundary  layer  which  grows 
from  500m  to  about  900m  in  4  hours,  and  has  a  Monin-Obukhov  length,  L,  growing 
from  100m  to  about  300m,  which  is  in  the  correct  range  to  obtain  longitudinal 
roll  convection.  We  estimate  that  critical  conditions  are  detectable  in  the 
rolls  for  about  1  -  1-1/2  hours,  and  also  that  average  statistics  need  to  be 
taken  over  a  period  of  about  1  hour  to  be  reliable;  these  experiments  force 
us  to  integrate  for  at  least  3  hours. 

Initial  conditions  for  the  model  runs  were  obtained  from  a 

one-dimensional,  steady-state  solution  using  the  artificial  cooling  of  Mason 
and  Sykes  (1982).  A  two-dimensional  perturbation  was  superimposed  to  initiate 
the  convective  rolls,  and  the  surface  temperature  was  fixed  at  the  initial 
value.  An  unfortunate  result  of  this  procedure  was  that  the  two  models, 
mixing  length  and  second-order  closure,  produce  different  steady-state 
profiles  from  the  same  external  parameters,  and  therefore  it  was  very 
difficult  to  obtain  a  very  precise  comparison  between  the  models.  The 
steady-state  initialization  is  important  because  imbalances  in  the  wind  field 
produce  inertial  oscillations  via  the  Coriolis  term  which  have  an  18  hour 
period,  and  would  consequently  affect  the  entire  run. 

The  initial  conditions  were  obtained  with  10  -  40m  for  the  mixing  length 
model,  and  A*  *  65m  for  the  second-order  closure  model.  Runs  were  made  using 
these  scales,  and  also  using  doubled  values,  i.e.,  80m  and  130m,  to 

investigate  sensitivity  to  closure.  All  runs  were  made  using  a  domain 
oriented  such  that  the  longitudinal  areas  of  the  rolls  lie  10°  to  the  left 
(facing  downwind)  of  the  geostrophic  flow  direction.  This  direction  was  found 
to  produce  the  maximum  momentum  transport  by  Mason  and  Sykes,  and  is 
consistent  with  field  observations. 

Figures  C.1  and  C.  2  show  a  sequence  of  realizations  of  the  flow  field 
from  the  two  models.  Vertical  velocity  component,  w,  and  temperature  contours 
are  shown.  The  S0C  model,  with  A*  »  65m,  produces  rolls  with  less  energy  than 
the  40m  mixing  length  model.  The  inversion  structure  and  internal  gravity 
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2-D  Roll  Calculations 

Run  Number  -  2102 
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FIGURE  C.l  :  FLOW  FIELDS  FROM  THE  TWO-DIMENSIONAL 
CLOSURE  MODEL  WITH  A*  =  65m.  HEAVY 
CONTOURS  ARE  VERTICAL  VELOCITY. 
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2-D  Roll  Calculations  -  Mixing  Length  Model 

Run  Humber  -  0011 
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FIGURE  C .2 :  FLOW  FIELDS  FROM  THE  MIXING  LENGTH 

MODEL  WITH  1  =  40m. 
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wave  field  is  much  better  resolved  in  SOC,  and  seems  to  show  a  transient 
generation  of  wave  trains  by  inversion  distortions  with  a  lifetime  of  about 
30  mins.  These  features  are  produced  by  a  large  amplitude  roll,  and  move 
along  the  inversion  with  the  roll.  A  wave  train  is  set  up,  propagating  energy 
upward,  until  the  roll  decays  and  the  source  of  energy  is  reduced.  The  wave 
train  then  detaches  from  the  boundary-layer,  and  moves  upward  with  the 
appropriate  group  velocity.. 

The  mixing  length  run  looks  different,  firstly  because  the  domain  is 
viewed  from  the  opposite  side,  so  that  the  direction  of  horizontal  propagation 
is  reversed,  and  secondly  because  the  rolls  are  more  intense,  producing 
smaller  scale  disturbances  in  the  stable  layer.  In  fact,  the  model  resolution 
around  the  inversion  appears  inadequate  for  this  run. 

Figure  C.3  shows  a  sequence  of  realizations  from  the  larger  parameterized 
scale  runs,  i.e.,  A*  -  130m,  and  a  mixing  length  of  80m.  The  mixing  length 
model  results  are  now  much  closer  to  the  A*  -  65m  pictures,  while  the  SOC 
results  show  very  much  weaker  rolls  and  wave  activity. 

The  time  and  space-average  heat  fluxes  from  both  SOC  runs  together  with 
the  one-dimensional  solutions  of  the  same  flows  are  shown  in  Figure  C .  ^ ,  and 
the  mean  temperature  profiles  in  Figure  C.5.  There  is  an  encouraging 

agreement  between  the  three  SOC  results,  although  the  partition  between 
resolved  and  parameterized  fluxes  changes  dramatically.  We  should  mention 
that  profiles  of  other  quantities  such  as  turbulence  energy  levels  or  momentum 
flow  would  look  very  different  for  the  various  models,  because  these 
quantities  contain  the  effects  of  internal  gravity  waves.  The  waves  would  be 
sensitive  to  domain  orientation  also,  as  shown  by  Mason  and  Sykes.  However, 

the  waves  do  not  transport  heat  flux,  and  therefore  the  entrainment  fluxes  can 

be  predicted  equal  by  the  different  models. 

A  similar  feature  can  be  seen  in  the  two  mixing  length  runs,  with 

resolved  flux  being  replaced  by  parameterized  flux  as  the  scale  is  increased 
(Figure  C.6).  The  magnitude  is  also  similar  to  the  SOC  results.  These 
results  cannot  be  readily  compared  directly  with  the  SOC  profiles;  because 
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2-D  Roll  Calculations 

Run  Humber  -  3003 


Time  *  10000.0  Sec 


Y  (  Km  ) 


Variable 

Units 

Increment 

Vertical  Velocity 

H/Sec 

0.25 

Temperature 

*C 

0.60 

FIGURE  C  . 3  :  CLOSURE  MODEL  WITH  A*  =  130m. 
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Model  Comparisons 
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FIGURE  C .4 ;  HEAT  FLUX  PROFILES  AVERAGED  FROM  t  =  9000s  to  t  =  12000s  FROM  THE 
DIFFERENT  MODELS. 
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FIGURE  C .5 :  AS  FIGURE  C.4,  BUT  FOR  TEMPERATURE  PROFILES 
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FIGURE  C . 6 :  HEAT  FLUX  PROFILES  FROM  t  =  9000s  TO  t  =  12000s  FROM  THE  TWO 
MIXING  LENGTH  RUNS. 
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the  surface  temperature  is  higher,  giving  a  larger  surface  heat  flux.  This  is 
the  problem  with  initial  conditions  which  was  referred  to  earlier. 

In  addition  to  providing  independent  calculations  of  the  turbulent  flow 
profiles,  the  two-dimensional  studies  also  provide  valuable  insight  into  the 
dynamical  mechanisms,  since  the  large  eddies  are  computed  explicitly.  We  can 
see  from  the  realizations  of  the  flow  field  from  the  SOC  model  that  there  are 
significant  vertical  displacements  of  the  inversion  produced  by  the  large 
eddies.  However,  the  inversion  does  not  suffer  any  large-scale  overturning  or 
breakdown,  so  there  is  no  entrainment  or  detrainment  due  to  large  eddies 
engulfing  stable  air  or  exchanging  mixed  layer  air  for  stable  air.  This  seems 
to  leave  two  mechanisms  for  entrainment;  small  scale  processes  at  the 
inversion  mixing  stable  air  downward,  or  large  eddies  providing  the  mixing  but 
in  some  way  not  involving  overturning.  The  simplest  mechanism  we  can 
postulate  is  that  the  large  eddy  carries  fluid  up  to  the  inversion,  where  it 
gains  some  heat  from  the  parameterized  mixing  processes,  then  the  fluid 

parcels  continue  in  the  circulation  and  are  brought  down  into  the 
boundary-layer.  Sample  fluid  particle  trajectories  from  the  80m  mixing  length 
run  are  shown  in  Figure  C.7;  particles  which  originate  near  the  surface  at 
the  base  of  an  updraft  and  are  carried  up  to  the  inversion,  are  invariably 
carried  back  down  again.  This  tends  to  confirm  the  simple  picture  of 
entrainment. 

There  is,  however,  a  problem  in  reconciling  this  simple  conceptual  model 
with  the  observation  that  the  entrainment  fluxes  can  be  almost'  entirely 
carried  by  the  large  eddies,  i.e.,  the  small-scale  mixing  is  not  necessary. 
We  must  therefore  seek  an  inviscid  entrainment  mechanism.  It  is  appropriate 
to  recall  here  that  we  are  dealing  with  a  boundary-layer  which  is  growing  due 
to  a  rising  mixed-layer  temperature  via  the  surface  heat  flux.  We  can 
consider  the  simple  "encroachment"  concept  of  Carson  and  Smith  (1971)),  which 
says  that  overlying  fluid  becomes  part  of  the  mixed  layer  as  the  temperature 
of  the  latter  rises  above  that  of  fluid  at  the  base  of  the  inversion.  This  is 
effectively  a  "no  motion"  concept,  and  does  not  predict  any  downward  heat 
flux;  stable  air  simply  becomes  part  of  the  mixed  layer.  This  cannot  explain 
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RUN  0013  PARTICLES 


FIGURE  C.7:  PARTICLE  TRAJECTORIES  FROM  THE  lc  =  80m  MIXING  LENGTH  RUN  FROM 
t  =  8000s  TO  t  =  12000s 
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our  observations,  and  so  we  suggest  a  model  which  takes  account  of  the 
dynamics  of  the  large  eddies.  If  we  consider  the  motion  of  these  eddies,  then 
there  is  an  inviscid  force  acting  at  the  inversion,  namely  the  dynamic 
pressure  gradient  force.  The  magnitude  of  this  force  will  be  related  to  w*/zi 
where  w#  is  the  convective  scaling  velocity;  which  characterizes  large  eddy 
velocities.  This  pressure  gradient  will  act  on  the  overlying  stable  air,  and 
if  we  consider  a  continuous  temperature  profile,  then  this  force  will  be 
capable  of  dominating  the  lowest  layer  whose  buoyancy  deficit  is  too  small  to 
resist  vertical  accelerations.  If  we  consider  a  layer  whose  temperature 
excess  over  the  mixed-layer  is  less  than  AG,  then  this  layer  will  be 
accelerated  by  the  pressure  gradient  and  can  be  brought  into  the  mixed-layer 
provided 
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The  downward  heat  flux  carried  by  large-scale  eddies  pulling  this  layer 
down  into  the  mixed-layer  will  be  roughly  w'AG,  where  w'  is  a  scaling  velocity 
for  the  vertical  velocities  near  the  inversion.  We  expect  w'  to  be  a  fraction 
of  w*  which  is  roughly  proportional  to  the  ratio  of  the  inversion  thickness  Ah 
to  z^  so  that  the  entrainment  flux  scales  as  • 


w*AG 


the  surface  heat  flux,  from  the  definition  of  w# .  We  therefore  have  a  simple 
mechanism  whereby  large  eddies  can  pull  fluid  from  the  base  of  the  inversion 
by  means  of  pressure  gradient  forces,  and  produce  plausible  estimates  of  mass 
and  heat  entrainment. 
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APPENDIX  D 


COMMENTS  ON  SCALAR  DIFFUSION  IN  THE  CONVECTIVE  BOUNDARY  LAYER* 


by 

W.  S.  Lewellen,  R.  I.  Sykes  and  S.  F.  Parker 

Wyngaard  (1984)  has  proposed  a  relatively  simple  parameterization  of 
diffusion  across  the  convective  boundary  layer  based  on  the  large  eddy 
simulation  (LES)  results  obtained  by  him  and  his  colleagues;  Wyngaard  and 
Brost  (1984),  Moeng  and  Wyngaard  (1984).  A  central  feature  of  this 
parameterization  is  the  incorporation  of  the  difference  between  the  effective 
eddy  diffusivity  for  scalar,  diffusion  down  from  the  top  of  the  convective 
layer  and  that  for  diffusion  up  from  the  bottom.  This  asymmetry  exhibited  by 
the  LES  results  for  the  scalar  transport  in  the  buoyantly  produced  turbulence 
driven  by  surface  heating  can  not  be  simulated  using  first-order  K  theory. 

In  this  comment  we  wish  to  address  two  questions  related  to  the 
asymmetric  diffusion.  First,  what  level  of  turbulence  closure  is  required  to 
exhibit  the  asymmetry  between  bottom-up  and  top-down  diffusivity?  Second,  how 
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sensitive  is  Wyngaard's  parameterization  to  the  details  of  the  asymmetric 
diffusion?  Not  surprisingly ,  we  find  that  a  level  of  closure  which  includes 
some  diffusion  of  the  second-order  scalar  correlations  will  produce  the 
asymmetry.  More  surprising,  we  find  that  Wyngaard's  parameterization  is  not 
sensitive  to  the  asymmetric  diffusion.  In  fact,  it  appears  that  precise 
symmetry  could  be  imposed  and  make  little  practical  difference  in  his  end 
results. 

Under  the  horizontally  homogeneous,  quasi-steady  conditions  assumed  for 
the  LES,  an  expression  for  K  in  buoyantly  produced  turbulence  may  be  written 
symbolically  as: 


^  \ 

/ 

K  -  -  w ' c '  /  (3C/3z)  -  Tl 

w f  w 1  - 
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where  the  time  scale  Ti  is  defined  equal  to 
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Under  the  same  conditions  the  potential  temperature- species  correlation 
appearing  in  Equation  D.1  may  be  written  as: 
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with  equal  to  the  ratio  of  c'0'  to  the  rate  of  dissipation  of  c'0'.  When 
Equations  D.1  and  D.2  are  combined 
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Equation  D.3  is  an  exact  expression  for  K  for  the  quasi-steady, 
horizontally  homogeneous,  buoyant  convective  layer  as  long  as  the  t's  remain 
exact.  If  the  t's  are  taken  as  properties  of  the  turbulence,  independent  of 
the  C  distribution  and  the  3rd  order  diffusion  terms  are  ignored  then  it  is 
clear  that  K  will  also  be  a  property  of  the  turbulence  only. 

Figures  8  and  10  of  Moeng  and  Wyngaard  (1984)  for  the  t's  as  given  by 

their  LES  results  suggest  that  for  this  problem  the  t's  may  be  taken  as 

properties  of  the  turbulence.  Thus,  the  asymmetry  appears  to  be  imposed  by 

the  turbulent  transport  of  either  c'0'  or  w'c'. 

The  answer  to  our  first  question  is  that  a  level  of  closure  that 
incorporates  turbulent  diffusion  of  the  second-order  species  correlations  is 
required  to  produce  a  proper  asymmetry  in  the  top- down- bottom- up  diffusion. 

Figure  D.1  compares  Wyngaard  and  Brost's  bottom-up  and  top-down 

dif f usivities  with  values  we  obtain  from  the  second-order  closure  model  of 
Lewellen  (1977)  which  uses  a  simple  gradient  diffusion  model  for  the  turbulent 
transport  of  both  c'0'  and  w'c'.  The  most  obvious  difference  between  .the  LES 
results  and  the  second-order  closure  results  is  the  singularity  which  appears 
in  near  z  =  .65z^  in  the  SOC  results,  but  doesn't  appear  in  the  LES 

results.  However,  this  is  apparently  more  a  result  of  uncertainties  in  the 
LES  result  than  it  is  a  true  difference  between  the  LES  and  the  SOC  results. 
Moeng  and  Wyngaard  (1984),  in  a  recalculation  of  the  LES,  present  results  for 
the  normalized  gradient  gb  which  go  thru  zero  a  little  above  z  =*  0.6z^.  Thus, 
their  more  recent  results  would  give  a  singularity  in  Kb  near  the  same  value 
of  z  as  given  by  the  SOC  results.  Since  the  turbulent  transport  terms  are 
divided  by  the  concentration  gradient  in  Equation  3,  relatively  small  errors 
in  the  turbulent  transport  can  lead  to  large  errors  in  K  in  the  middle  of  the 
convective  layer  where  the  gradient  is  very  small.  We  conclude  that  the 
asymmetry  between  Kb  and  Kt  exhibited  by  _the  SOC  model  is  qualitatively 
similar  to  that  of  the  LES.  Presumably  we  could  use  terms  of  the  type 
proposed  by  Zeman  and  Lumley  (1976)  to  adjust  the  modeled  turbulent  transport 
terms  and  make  our  profiles  of  Kt  and  come  closer  to  the  LES.  As  more 
detailed  LES  results  become  available  this  may  be  desirable. 
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Dimensionless  Diffusivities 


Kt.  Kb 


FIGURE  D.l:  COMPARISON  OF  THE  DIMENSIONLESS  SCALAR  DIFFUSIVITIES  AS  OBTAINED 
BY  THE  LES  OF  WYNGAARD  AND  BROST  (1984)  WITH  THOSE  FROM  A  SIMPLE 
SECOND-ORDER  CLOSURE  MODEL. 
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The  next  question  we  address  is  how  important  is  this  K  asymmetry  in 
Wyngaard's  proposed  parameterization.  We  approach  this  question  by  repeating 
Wyngaard's  analysis  with  the  top-down  dimensionless  gradient,  gt,  multiplied 
by  a  factor  K-j  .  The  parameterizations  provided  in  Equations  21  and  23  are 
then  modified  by  replacing  cw-j  and  we  with  cw^/K^  and  we/K^  respectively.  The 
only  other  change  is  Equation  29  for  At.  It  becomes 


(D.  4) 


These  changes  force  two  types  of  changes  in  the  resulting  unmixed  layer 
profile.  First,  h^  is  moved  a  little  further  from  z^  as  is  decreased,  and 
second  the  gradient  of  C  in  the  interfacial  region  is  increased  (or  decreased) 
as  K-j  is  increased  (or  decreased).  To  make  this  sensitivity  more  specific, 
consider  Wyngaard's  boundary  conditions  for  humidity-like  scalars.  Then  his 
Equations  26  to  31  yield 


(D.5) 


If  the  top-down  and  bottom-up  diffusivities  were  made  symmetric  by  making 
K.j  =■  0.4  then  cw^  would  be  decreased  by  *  20%  for  typical  PBL  conditions  by 
virtue  of  WgA^/w*  increasing  from  *  0.2  to  «  0.4.  However,  since  we  cannot  in 
general  be  estimated  to  within  an  accuracy  of  20%  this  change  is  of  little 
practical  significance  until  better  parameterizations  of  the  entrainment 
across  the  inversion  are  available. 

Wyngaard's  section  4  is  devoted  to  arguing  that  his  closure  is  valid  for 
a  wide  variety  of  typical  conditions.  We  believe  that  this  is  true,  but  not 
because  the  LES  results  for  gt  are  valid  for  a'll  these  conditions.  It  is  true 
because  the  critical  features  which  allow  the  parameterization  to  work  are 
that  K  go  smoothly  to  a  reasonable  representation  of  the  surface  layer 
dynamics  in  the  lower  part  of  the  layer  and  that  K  go  smoothly  to  a  small 
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value  within  the  capping  layer.  We  believe  the  LES  results  for  gt  in  the 
upper  part  of  the  boundary  layer  are  likely  to  be  sensitive  to  such  things  as 
strength  of  the  inversion  and  non-linearity  in  the  flux  profile,  but  this 
sensitivity  is  of  little  consequence  until  a  much  better  parameterization  of 
the  inversion  layer  is  available. 

In  Summary,  the  asymmetry  in  scalar  diffusion  within  the  quasi-steady, 
homogeneous,  convective  layer  is  an  interesting  feature  which  can  be  used  to 
test  diffusion  models.  Second-order  closure  has  a  definite  advantage  over 
lower  level  closure  in  simulating  this  feature.  However,  further  improvements 
in  the  parameterization  of  entrainment  across  the  inversion  are  necessary 
before  accurate  simulation  of  this  interesting  feature  is  very  significant  in 
practical  mixed-layer  parameterizations  such  as  that  proposed  by  Wyngaard. 
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